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This report summarizes the work carried out in the past three 
years under the Lunar Surface Electric Properties (SEP) Project. 

It concerns primarily with the theory of rad io- frequency interfero- 
metry used in geophysical subsurface probing. Parts of this report 
have been published in the following articles: 

1. Kong, "Electromagnetic Fields Due to Dipole Antennas over 
Stratified Anisotropic Media," Geophysics, Vol. 37, pp. 985- 
996, 1972. 

2. Shen, Tsang, and Kong, "Multifrequency Excitation of a Wire 
Antenna for an Invariant Radiation Pattern," IEEE Trans, on 
Ant. and Prop., 1972. 

3. Kong, Tsang, and Simmons, "Lunar Subsurface Probing with Radio 
Frequency Interferometry," URSI Symposium, 1972. 

4. Tsang, Kong, and Simmons, "Interference Patterns of a Hori- 
zontal Electric Dipole over Layered Dielectric Media," 

J. Geophysics Research, 1973. 

5. Kong, Tsang, and Simmons, "Geophysical Subsurface Probing with 
Radio Frequency I nt er f er o me tr y, ** IEEE Trans, on. Ant. and Prop., 
July, 1974. 

6. Tsang, Brown, Kong, and Simmons, "Numerical Evaluation of Electro 
magnetic Fields due to Dipole Antennas in the Presence of 
Stratified Media,” J. Geophysics Research, 1974. 

Contributions are also made by L. Tsang, John Mallick, Winston Chan, 
Paul Palmer through their S.B. Theses, and by L. Tsang through his 


S.M. Thesis. 
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Chapter 1 . 
INTRODUCTION 


The subsurface of planetary bodies, including the earth and 
its moon, can be examined with electromagnetic waves. In 1955, 
the radio-frequency interference fringes method was used on geo- 
physical prospection of underground water in Eqyptian desert 
(El Said, 1955a, 1956b). This report concerns with the experi- 
ment developed for Apollo 17 mission to measure the subsurface 
electromagnetic properties of the moon (Simmons et al, 1973). 

In the experiment, a transmitting antenna consisting of a pair 
of orthogonal dipoles is laid directly on the lunar surface. The 
antenna radiates sequentially at frequencies of 1, 2.1, 4, 8.1, 

16 and 32.1 MHz in a time window of 100 millisecond for each fre- 
quency. The time window is sufficiently long that the experiment 
is effectively a continuous wave experiment. A receiving antenna 
consisting of three orthogonal loops is mounted on the lunar roving 
vehicle. As the LRV traverses the lunar surface, the strengths 
of magnetic field components are measured as a function of distance 
from the transmitting antenna. The data are recorded on magnetic 
tapes and returned to earth for analysis. The interference pat- 
terns of the field-distance plot contain information about the 
electromagnetic properties of the lunar subsurface. To test both 
equipment and theory, experiments have been performed on three 
glaciers (Strangway et al 1974) in Switzerland, Canada, and Alaska 
and on a scale model tank (Rossiter et al, 1974) in laboratory. 
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In this report, we present the theoretical basis of the 
experiment. The mathematical model is a stratified n-layer me- 
dium. Each layer is bounded by plane boundaries and possesses 
different electric permittivity, magnetic permeability, and 
thickness. Although dipole radiation in the presence of stra- 
tified media has been studied extensively (Sommerfeld 1949, 
Brekhovikih 1960, Banos 1962, Gudmaudsen 1972, Wait 1970, Felsen 
and Marcuvitz 1973, Kong 1972, Annan 1973), a proper account for 
the interference fringes method is still lacking. The geometri- 
cal optics approximation was the first approach used in calcu- 
lating interference patterns (El Said 1956, Annan 1973, Tsang 
et al 1973). With the use of the reflection coefficient formu- 
lation (Kong 1972) , all field components can be expressed in 
integral forms with a single variable of integration. In order 
to obtain explicit expressions for the six field components, the 
following three different approaches are used to evaluate the 
integrals: 1) geometrical optics approximation, 2) modal ap- 

proach, and 3) direct numerical integration. Scattering effects 
are Also discussed: the advantages and disadvantages of the’ 

various methods and their respective regions of validity are 
compared in the last Chapter. All numerical results for aniso- 
tropic as well as isotropic cases are presented in Appendices. 
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THEORETICAL FORMULATION 

2.1 INTRODUCTION 

The problem of radiation of a dipole source in the presence 
of stratified media has been extensively studied with application 
to geophysical explorations . An excellent review on the half-space 
case is contained in the book by Sommerfeld (1949) and in the- 
monograph by Banos (19 66) . Propagation and Radiation in stratified 
media are treated by. Wait (1970) amd Ward (1967) . Wolf (1946) , and 
Bhattacharya (1963) considered the case of dipoles on two layer 
earth. Wait (1951, 1953) solved the problem of electrical and 
magnetic dipoles over a stratified isotropic medium. The case of 
an anisotropic half-space was studied by Chetaev (1963) and Wait 
(1966a). Praus (1965), Sinha and Bhattacharya (1967) and Sinha 
(1968, 1969) treated electric and magnetic dipoles over a two layer 
anisotropic earth. Wait (1966b) formally solved the case of a 
horizontal dipole over a stratified anisotropic medium. All these 
works are carried out by means of Sommerfeld* s Hertzian potential 
functions, and the primary interest is concentrated in the limits 
of high conductivity. Magnetic properties are almost entirely 
neglected? mainly because such model studies assume principal 
applications to the earth where the permeability is nearly equal 
to that of vacuum and the electric conductivity dominates at low 
frequencies. In events of other celestial bodies, such as the 
moon where the lack of moisture renders very low conductivity to 
the medium, a study of contributions due to all electric and 
magnetic properties then becomes important. 
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This Chapter is devoted to the case of radiation of various 
dipole sources in the presence of a stratified anisotropic media. 

The anisotropic medium is uniaxial and possesses both tensor 
permittivities and permeabilities. The principal axes are all 
perpendicular to the boundaries separating different media. 

Solutions to the problem are facilitated by decomposing a general 
wave field into TM and TE modes, employing the concept of propaga- 
tion matrices and expressing the reflection coefficients in terms 
of continuous fractions . The primary excitation is separated entirely 
from contributions due to the medium. The reflection coefficients 
depend solely on the geometrical configurations as well as the 
physical properties of the stratified medium. 

In studying the theory of electromagnetic wave propagation, 
it has been appreciated both classically and quantum mechanically 
(Kong, 1970) that introduction of a potential function is not 
necessary and sometimes complicates the algebra especially when 
anisotropic media are involved. With recognition of the fact 
that outside any source, two scalar functions are sufficient to 
determine all field quantities, two components of the field 
vectors can be chosen as the fundamental scalar functions. In 
our case, the preferred field components for TM and TE decomposition 
are clearly those along the principal axis and to boundaries of 
stratification. With the aid of propagation matrices (Kong, 1971), 
wave amplitudes in any region are easily calculated in terms of 
those in any other region. Writing in the form of continuous 
fractions, we obtain a closed form solution for the reflection 
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coefficients. All field components are expressed in terms of 
integrals which are ready for direct numerical evaluation. A 
discussion is given for the various special cases. 
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2.2 TRANSVERSE ELECTRIC AND MAGNETIC WAVES 

Governing equations for electromagnetic fields in a region 
outside any source are the Maxwell's source free equations. 


Vx E — iuup • H (2.1a) 

V x H = - iuui . E (2.1b) 


where in (2.1a), y is the permeability tensor of the media. The 
tensor e in (2.1b) contains information about the dielectric 
constant and the conductivity of the medium. e = e' + ie" , 
where e 1 is the permittivity tensor, and 'e" is related to the 
conductivity tensor a by e“ = a/w. Time harmonic excitations 
with time dependence exp(-iwt) have been assumed. The tensors 
c and y can be represented by hermitian matrices. In our case, 
we consider media which are uniaxially anisotropic, where 




and 


f M 




u 





(2.2a) 


(2.2b) 
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We employ cylindrical coordinates, and the plane transverse to 

A 

the 2 axis is characterized by p and <p . Longitudinal elec- 
tric and magnetic components E and H are used to derive TE 
and TM waves. The wave equations to be satisfied by E and H, 
are immediately derived from equations {2.1) and (2.2). If we 
take the z component of (2.1b) in view of e given by (2.2a), 
employ (2.1a) to eliminate transverse magnetic field components, 

and use the fact that V • E = (1 - a)3E /3z, the equation for 

z 

E z is 


(V t + a 2 + ^ a ) E z “ 0 


By the same token, we obtain the wave equation for H, 


(V/ 2 + b -~r + k 2 b) H = 0 

t _ 2 Z 

Sz 


In equations (2.3a) and (2.3b), 


k = to/ 


a = e / e 

Z 


b = ^ / u 


iil.iij.i-i- 

Pa ” p 2 a„ 2 
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is the transverse Laplacian operator expressed in cylindrical 
coordinates. It is seen from Eq.(3) that E and H are 

Z Z 

decoupled, which would not be true if the e abd y tensors 
possess off-diagonal elements. A unique decomposition of the 
total wave into a transverse magnetic field (TM) mode derivable 
from E and a transverse electric field (TE) mode derivable 

Z 

from H is therefore plausible. We note that a pair of vector 
z 

wave equations can be derived from (2.1) and (2.2): 

V 2 E + k 2 E + (a - l)k 2 E z + (a - 1)V(3E /3z) » 0 (2.7a) 

z z 

V 2 H + k 2 H + (b - l)k 2 H z + (b - 1)V(3H /3z) = 0 (2.7b) 

z z 


Eq. (2.7a) is the vector wave equation for the electric fields of 

TM waves where H =0, and (2.7b) is the wave equation for the 

z 

A 

magnetic fields of TE waves, where E = 0. The z component of 

z 

the two vector equations (2.7) gives rise to equation (2.3). 

Solutions of E and H to the wave eouation (2.3) in 

z z 

cylindrical coordinates are well known. As a consequence of the 

Maxwell equations (2.1), all transverse-electric- and magnetic- 

field components can be expressed in terms of the longitudinal 

components E and H which, respectively, characterize the 
z z 

TM and the TE waves. In our problems we are interested in wave 

A 

solutions which are outgoing in p direction and traveling or 

A 

standing in z direction. Therefore, we obtain, for a fixed 
separation constant n , 
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ak ( ^ -ik^z ik^z , _ M 

i £ ( _ Ae 2 + Be ) H n (1) (k p p) s n (cp> 

P 


, (€> ... (e) (e) 

ak — ik z lk z v tm* 

i_|_(. Ae * + Be Z ) H n (1) (k p p)S n ™ ^ 

V 


ik (e, z ik (e, z 

CA(k Z + B(k ). Z ] H (k p) S ™ (cp) 

p p n p 


ik^z ik^z , 

- i - a ~~ ( Ae 2 + Be 2 ) hJ (* p P>S n (cp) 

% P 


-ik^z ik^z , TM 

i ^ ( Ae 2 + Be 2 > H n 1 (k p P )S n (cp) 

P 


0 
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-TE 

E 


co 


dk 


-co 


-ik (m) z ik (m) z 

i^(Ca 2 + Be 2 > H' > l >(k 0)S n TE '( V ) 

- z n p n 


-l 


..(m) .,(ra) 

hiDii z 2 ^ z 2 (1) 1 TE 

( Ce + De Z ) H^ 1 ' (k p p)S n TE (cp) 


(2,8c) 


i Sl 


-te 


OD 

f dk 
J— P 


-oo 


(m) (m) 

z Z lk z 2 (1) ' TE 

k ( -Ce Z + De ) H n 1 (k p) <<P> 

r 


bk 


(in) 


— -ik (m) Z ik ^ )z (1) TE 1 

.2 ( -Ce 11C z 2 + De 2 ) H U, (k tfS™ ( cp) 

k p n p xx 


— ik z ik z 

1,K z 2 2 z m TE i 

[C(k ) e + D (k )e ] H ' (k p) S_ (cp) I 

V P p n P n y 


(2.8c 



2-9 


where superscripts TM and TE denote, respectively, TM and TE waves. 
We note that if the integrands for E z and H z are denoted, 
respectively, by E z (k p ) and H z (k p ) such that 


E 


TM 

Z 



dk E (k ) 
p z p 


and 




dk P H z ( y ' 


then the integrands of the transverse components are related to 
E (kp) and H <kp) by the following relations! 

2 Z 


w 


TM 



9 W /3z ' 


w 


TM 


= -l 


awe 


^ x E„ (k ) , 


and 


(2.9a) 


W* - ; *t| 


t p 


dH (k )/3z 
Z P 


, E t (k p ) TE = * H z (k p ) 

P 


(2.9b) 


where 


A _ A A a a1 

E = zE„ , H - zH , and V. ■ *» p •*— . A— ■&- 
z z ' z z t M ^p + ¥ p 3p 


The fact that TM waves are extraordinary waves in the medium is 

signified in equation (2.8) by the superscript (e) on the z- 

( 0 ) ■ * 

directed propagation constant k z , which satisfies the disper 
sion relation 

k (e) - <k 2 - k 2 /a) 1/2 . 
z M 


(2.10a) 
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TE waves are derived from H and satisfy the dispersion relation 

Z 

k (m) = Oc 2 - k 2 /b) 1/2 , (2.10b) 

2 P 

where the superscript (m) indicates the effect of magnetic 
anisotropy. In (2.8) , the first element of the column matrices 

A A 

denotes p component , the second element of <}> component, and 
the third element the z-conponent. The Hankel functions 
of the first kind and nth order represent outgoing waves in 
direction due to our choice of the time dependence exp(— iwt) • 

S (4>) stands for sinusoidal functions of . Primes on H ^ (k p) 
n n p 

and S {$) denote differentiation with respect to the arguments, 
n 

The k dependent functions A, B, C and D are to be determined 
P 

by the appropriate boundary conditions. 

2.3 PRIMARY EXCITATION 

The explicit solution to the problem of dipole radiation 
over stratified medium depends on field excitations of the source, 
and the geometrical configuration and physical constituents of the 
medium. In the absence of the stratified medium, the solution of 
electromagnetic fields in an isotropic medium due to a dipole 
antenna, which we refer to as the primary excitation, is well-known 
(Adler, Chu and Fano, 1960). The solution is usually written in 
chperical coordinates. It can be transformed into cylindrical 
coordinates and represented by Hankel functions in the integral 
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form. Writing in the general form, we have 


E 


co e 

f dk E (k ){ 

J _ DO p ] -1 


ik z a 
z ' 


z »> P 

-oo K 


p ) - iJt z z 
e “ 


H (1) (k p) S ™ (cp) 
n p n ^ 


z > 0 (2.11) 

z < 0 


oo 

H = r dk H (k ) 
2 ioo p ° p 


f xk z \ 
' z 1 
e 


-ik z 
z 


e 


H 


( 1 ) 


TE 


n 


<V> V <*> z < 0 


z > 0 (2.12) 


where E and H characterize the structure and excitation of 
o o 

the dipole. All field components follow from Eq. (2.8) 

with B=D=0, A = E , C=H for z < 0, and A = C = 0, 

o o 

B = E , D = H for z a 0. 
o o 


For the elementary dipoles under consideration, we obtain 

for: 

TM 

1) Vertical electric dipole: n = 0, (cp) = 1 


E 

o 


ii* p 3 

8TTojek 

z 


(2.13a) 


H 


o 


0 


(2.13b) 
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2) Horizontal electric dipole along x-direction 
Ilk ^ 

P > TM 

E = + i — z < 0 S, = coscp 

o — 8TTiue 1 


(2.14a) 


Ilk, 


H o 1 8 n k 


_ TE 

S, = - smcp 


(2.14b) 


TE 

3) Vertical magnetic dipole: n = 0, (cp) = 1 


IAk 


H = 
O 


- x 


8TTk 


(2.15a) 


-JE = 0 
o 


(2.15b) 


4) Horizontal magnetic dipole along x-direction 

2 


IAk, 


H = + — — 
o 8tt 


z < 0 


TM 


= COS cp 


(2.16a) 


IAiupk ( 
E o " 8 rr k _ 


0 TE 

= - sxncp 


(2,16b) 


In the following, we derive Eq. (2.14) from well-known potential 
solutions for the dipole. Similar derivation, comparatively sim- 
pler, applies to Eqs. (2.13), (2.15) and (2.16). The vector poten- 

tial solution for the horizontal electric dipole is 
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_ A A 

A = A(pcos<j> - <J>sin<f>) 
where A = e^^ r /r 


The electromagnetic fields are obtained from 
H = V x A 

and 

E « i ig {V(V • A) + k 2 A } 

Using the identity (Sommerfeld, 1949) 


e ikz _ i f 
~ 7 J 


dk 


k 

P 


H 


( 1 ) 


(k pP ) 


+ik z 
- z 


the field components can all be written in the integral form. 
For the z-component. 


H »= - sincfi 5 A/dp 
z 

1 2 

E c i — r ccs<f> 3 A/3p3 z 
Z u)£ 


The results are Eq. (2.14). 

In Eqs. (2,13) - (2.16) I is the current that drives the di- 
pole, 1 the equivalent lingth of the electric dipole, A the area 
of the current loop that constitutes the magnetric dipole. Hori- 
zontal dipoles can be simply obtained by a rotation of coordinates, 
which ammounts to change cos<t> to sin<|> and sin4> to - cos4> . 

We note that a vertical electric dipole excites TM wave only and 
a vertical magnetic dipole excites TE wave only, both involve 



2-14 


Hankel functions of zeroth order; whereas horizontal dipoles excite 
both TM and TE waves and require Hankel functions of first 
order. An arbitrary oriented dipole can be treated as a linear 

A a a 

combination of three dipoles along x, y, and z axes which are 
just described. 

2.4 DIPOLE ANTENNAS OVER STRATIFIED ANISOTROPIC MEDIA 


Giometrical configuration of the problem is shown in Fig. 2.1. 
There are n slab regions, and the last region is numbered t 
instead of n + 1, for the sake of simplifying writings. In each 
region labelled i, solutions of electromagnetic field components 
take the form Eq. (2.8) with all quantities subscripted by i. In 
the Oth region where we have the antennas, Aq * and C^ * Hg 
which are known from (2.13) - (2.16) for the three types of anten- 
nas under consideration. In the last region namely region, t, it 
is semi-infinite and we do not expect reflected waves, therefore, 

B t = D t = 0. 

Boundary conditions at all interfaces are that all tagential 
electromagnetic field components must be continuous for all p and 
<f>. Consider the boundary at z *= -d^, the continuity of tangen- 
tial electric fields and the continuity of tangential magnetic 
fields yields, for the TM waves, 


- ji, (e)j 

{" A,e 1 +B. e 


-ik. (e) d. 


xz 


X > 


ik <3 

v / _ a « (i+1) z i ~ JV (i+3W“3 

K (i+l)z t A i+l e + e 11+1,2 i 


(e) 


te) 


ik, ? ' d, 

„ (i+l)z t x 

i+1 e > 


(2.17a) 
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ik 1 , ,e) d ) -ik. (e, d. 

e i « A i e + B ± « h 


ilc ^ _j% (®) j 

r A e < i+1)z d± + B e ik (i+l)z d i. 

e i+l( i+I e + B i+1 e > 


(2.17b) 


and for the TE waves 


ik. (m) d. -ik. (n0 d. 

y i { C i e + D i e 


„ rr 1K (i+l)z d i 

^i+l { C i+1 e + D i+l e 


i t (m) ,, 

lk (i+l)z d i 


> 


(2.18a) 


ik. (m, d. {m >d, 

k iz ^ ^i e + D ie 1Z i) 


ik <m) d -it (m> h 

-C.„ e < i+1 >* 1 + D i+1 e < r+1 >* S 


(i+1) z 1 i+1 


(2.18b) 


The reason that we can treat the TM and TE cases separately is 
because that 1) for vertical dipole case, only TM or TE is 
excited; 2) for horizontal dipole case, although the total tan- 
gential field components consist of both TM . and TE waves, the 
coefficients of and H^ 1 ' separates TM and TE cases. 

We now illustrate the derivation for the upward propagation 
matrix M^ + 1 for TM cases. From (2.17a) and (2.17b) it is 
straight-forward to solve for A^exp ik^^d^ and B^exp - ik^^d^ 
in terms of + 1 ^md B i + 1* De fi ne 

(e) 

= A. exp (ik. d. ) 

1 * 1Z 1 


a. 

i 


(2.19a) 
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b. - B . exp ( - ikf e) d. ) 
X x c XZ X ' 


We have 


= \ {e <+> ^ +1 e (+) i < e , a i+1 + C (-)| +1 e(-). (e, b. +1 ) 


b i - f <e(->J +1 e ( + )i <e ’ a i +1 + E(+)i +1 e <->i (e,b i+i> 


The upward propagation matrix for TM waves from (i + l)th 
to ith region is defined to be 


e {+) i+ j 1 e (+) 


e (-) 1+1 e (-) (e) 
1 1 


« i+1 

x 


1 

2 


ef -) 1 * 1 e (+) ( ®> 


e (+) 1 ^ 1 e (-) 


J 


where 


€ <±> q 


P _ 


= ( 


.(e) 


e 

e - . (e) 
q k 

qz 


e( i )( ! > ‘ exp Ct ik uii)z <a i ’ d i + i>l 


(2. 19b) 

(2.20a) 

(2.20b) 

region 

( 2 . 21 ) 
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(e) = r 

iz 


k. - k /a 
x o' 


Equation (2.20) can be written as 


fa. ) 
1 

= M i+1 

a i+l 


1 

• 

b 


b. , 

i- J 


L 1+1 , 


( 2 . 22 ) 


TM TM 

Note that, with R denoting reflection coefficients and T 

denoting transmission coefficients, we can write 


a 0 = E 0 exp (ik< z’ d c 


(2.23a) 


b 0 = R ™ E 0 eXp ( " ik< z* d 0 ) (2.23b) 

a t = T™ E Q exp (ik ( ® } d fc ) (2.23c) 

and that b fc = 0 because there is no reflected wave in the last 
region. The parameter d^ in (2.23c) is introduced for convenience, 
it does not correspond to any distance and is always multiplied by 
e(-)^ to yield exp(ik^ e ^ d ). The definition of the propagation 
matrix is a useful one. Once wave amplitudes in any region are 
known, those in the regions above this one are all determined by 
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(2.22). Thus the propagation matrix (2.21) propagates wave 
amplitudes upward. We can also define a propagation matrix M 
which propagates wave amplitudes downward. 


M i+1 


1 

2 


, . i , . (e) 

e(+) i+ i e(-) i 


e(-) ^i e (+) 


(e) 

i 





r " 




a i+l 

■ M i« 

r a i 

and 


b i 

* 



It is easily shown that 


« i t 1 «i+i= 1 


Following a parallel analysis of the above, we define a 
downward propagation matrix for TE waves . 


N.\ 

i+l 


1 

2 


. . i . x (m) 
U( + ) i+ l e(-) . 


e(+) ( i’ 




^ (+) i+l 6(+) 4 


(m) 

i 


i 

i + 1 


(2.24a) 


(2.24b) 


(2.25) 


(2.26a) 
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e i+l. ] 

= N. 1 

l+l 

c. 1 
1 

a 1 

A 

d i + l . 

) 

Cl . 

1 J 


(2.26b) 


i+1 


An upward propagation matrix N ^ is defined as the inverse of 


N iii 


r p(+) i t 1 e (+) <”> 


M 1+ l 1 
N i - 2 




1+1 e (+) (m) 

V(-) i e(+) i 


U(-) i J 1 e(-> 


]i(+) i l 1 e(-) ( J° 


(2.27a) 





c i 

= N 1 * 1 

c i+l 



k d i+l , 


(2. 27t ) 


In (2.26) and (2.27) 


M<±>$ 


v (*) 

Hi - k (m) 

qz 


(2.28) 


c. = 


C. exp(ikf m ^ d. ) 

i r 1Z i 


(2.29a) 


d. = 


D. exp( - ik d. ) 

1 1Z 1 


(2.29b) 


(m) 

iz 


- -/k. 2 - 

1 


k 2 /b 


(2.29c) 
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Also c Q = H Q , d Q » R TE H q , c t = T^Hq and d t = 0. We note 
that for vertical magnetic or electric dipoles, only the TM waves 
or the TE waves, respectively, are excited. In the case of either 
a horizontal electric dipole or a horizontal magnetic dipole, both 
TM and TE waves are excited. d 

2.5 REFLECTION COEFFICIENTS 

In the interferometry method, our primary interest is the re- 
flected wave fields. From the preceding section, we have established 
that when the wave amplitudes in region 0 are known, solutions in 
any other region can be determined by using the downward propagation 
matrices (2.24) and (2.26). In this section, we derive for the re- 
flection coefficients, a formula which is expressed in continuous 
fractions. We observe from Eq. (2.23) that b fc /a t = 0 and b Q /aQ 
gives rise to the reflection coefficients R . And from (2.22) , 
an expression for b^/a^ terms of b^ + ^/a^ + i can ^* e eas Hy 

extablished. In view of (2.21), (2.22) gives 



e (+) 
e(-) 


i+1 

i 

i+1 

i 


( 1 


/ \ i+1 / / \ i+l , \ i+1 / , \ i+1 

e (+) i / e(-) i - e(-) ± / e(+) ± 

•M 1 ? / s(-) 1 * 1 + (e(-) ( *>/(b i+1 / a . tl ) 


(2.30) 

Making use of (2.30), we obtain a formula in continuous fraction 

TM 

for the reflection coefficient R 
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R™ = exp(i2k (e) d„ ) ^ 

" 2 ° €<-)i 


(1 - 


«(*)p / e <->0 - / e <+>0 | 

I e(+)i / c(-)i 


+■ exp(i2k ( ®’ 


d 0» 


e(+) 
e (-) 


[1 - 


e{+)l A e(-)* - e (-)^ / e (+) ^ ; 

I e(+)J / e<-)* 


+ exp(i2k (e) (d 
* nz n 


e (-) 


Vi> > 


n 


®< + >n 


,]}. ( 2 . 31 ) 


Likewise, we obtain the reflection coefficient for TE waves 


_TE # . (m) j % 

R = exp(r2k a ) 
n z o 


n(+) 

u(-); 


{ i - 


H<+>0 / n(-)J - u(-)q / U(+)g 


r 


u( + )J / u<->0 


+ exp(i2k^™^ (dj 


M(+>i 

d o» — 72 C 1 
u(-) n 


P(+)^/n(-)i - u(-)^/u<+)J j 

[ " ~2~ 2 
y(+)l / u(“>i 


(m) 


U(-) 


+ .... + exp(12k nz (d n - d^)) 


n 


^< + >n 


]} 


( 2 . 32 ) 


mw rpp 

The subscript n on R and R denotes the number of layers 


involved 
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2 . 6 SUMMARY 


With the reflection coefficients determined in the last 
section, we can now summarize the formulas for all field quantities 
in region 0, where interference patterns are calculated. We have 
decomposed total wave fields into a summation of the TM and TE 
wave modes. 


E = 


-TM 

E 


-TE 
+ E 


(2.33) 


-TM 
H = H 


-TE 
+ H 


(2.34) 


The TM and TE solutions are 
1) For a vertical electric dipole 


-TM 

E 


= J dk (- — — ) 

J p 8 TTtue 


_ (*> (e) 

„ +ik z lk z .... 

ik (+ e 2 + R E 2 ) H. (k p) 

p - Op 


-oo 


. . je) 
+ik z 




2 / ” z . «tm _ z\ „ (l) 


(e) 


( e 


+ e * ) H '*'(k P) 

u p 




J 


-TM 

H 


-co 


II 

p ' 8 tt 


dk. (- tf) 

% 


(2.35a) 


k 2 +i^'z „„ ik W, z .... 

ipj ( « z + R™ e Z ) H 0 (1) <k p) 


(2.35b) 
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e te 



2) For a vertical magnetic dipole 


-TE 

E. 


oo 

■J 

-OO 




2 (**>l (vt'l 

k . +ik 2 ^ ikz 

. p , — z _TE z X „ (1) /I \ 
rs*) ( e + R e ) H q (k p) 

z 


J 

(2.36a) 


/' _ +ik z „„ ik z . , 

• i 2 , . — z _TE z . (1) . . 

ik (+e + R e ) H (k p) 

P Op 


aTE - L dk P 


-CO 


v 3 +i^ w z iiK 

k P , ilk z Z TE z _ (1) . 

— (e + R e ) H q (k p) 

rw ^ 


J 


-TM -TM 
E = H =0 


(2.36b) 


3) For a horizontal electric dipole along x-direction 


i 



2-74 


E™ 


oo 

■ L d \ (i 


. II 
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p 8rruoe 


(&) Ce) 

. +ik z _ M ik z . , ^ 

iak ® k (e z - r™ e 2 ) (k p)cosJ 


z p 
(e) 


k' w +ik~'z m> , ik Ce z ... 

-a — — (e 2 - R T e 2 )H. (k p)sincp 

p 1 P 


\ P 


2 ±^2 Z TM ^Z 2 (1) 

k 2 (+e - R™ e '(k p)coscp 


J 


(2.37a) 


H™ 


oo _ , 

■ I •% <-§^> 


-oo 


+ik (e Z ik fe) 2 

a z TM z. „ (1),. . . 

- — (+e - R e )H^ (k^pjsincp 


+ik^ z ilc^z 

-ak {+e 2 - R T e 2 )H (k p)coscp 

0 l o 


J 


(2.37b) 


-TE 
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oo 

-J 

-co 


dk 

p 8 tt 


(«) (n>) 

k +ik z __ ik z . 

^ (• Z + »**• > H 1 ( ’‘kp^oo., 

z 


kp +i^ TE ik ^ Z (1)’ 

b^(e +R e )H X ( (k p p)sincp 

2 




(2.37c) 



2-25 


-oo 


f Dc^z - . , 

ibk (+e~" Z + R TE e Z )H_ ( (k p) sin<p 
D X 0 


S TB f°° II . 

H = j (1 8^ > 


j. , M 

+ik z 


ik W z 


— (+e z + R TE e z )H- *■*■* (k pjcosq, 

0 “ 1 P 


k^ ilk z ik z . 

(e Z + R TE e Z jH^^k p) sincp 
K z p 


J 


(2.37a) 


4) For a horizontal magnetic dipole along x-direction 


oo 

=™ = J <- V 


“OO 


..,<*> . «) 

+ik z ik z . , 

iak (+e 2 + R TE e Z )H ± ^ } (k p) sincp 


j. ,ce) . ,Ce) 

. +ik z m _ xk z . 

ia, l “ z _TE z » (1) .. \ 
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2 . < c> te.) 
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(2.38a) 


,T» 


f» 

J _ p 8 tt 


-oo 


, <e) ce.) 

_ +ik z m _, ik z x 

- + R e ) H- (k p)coscp 

V 1 p 


ak j-ilc e, z ilc*z nY> 

i_ ^(e Z + R e Z ^ (k p)sincp 
z p 


(2.38b) 
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b 

P 


(±e 


+xk z 
— z 


+ R 


TE 


. 

ik z 
e 2 ) H. 
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(k p) cos 
P 




-TE 

E 


I dk 

-m 
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8tt 
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±ik z TE i* 2 (1) 

(+e 2 + R TE e 2 )H 1 (1) 
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-CO 


i b K<"> k ,-e— +RTEfi <V> C ° S * 
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<m> M 


bk'" 1 ' ±i* 2 TE i\z (1) 

— ^ — (-e + R TE e )H 1 (k p) sincp 


-X 
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0 +ik z ik z 

1 2 tJ—~ 2 . « TE z 

k (+e + R e 
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(2 . 38d) 



2-27 


2.7 DISCUSSIONS 


The problem of radiation of various dipole antennas over a 
stratified anisotropic medium has been solved. In view of the 
general formalism as presented in sections V and VI, we can make 
the following observations: 

1) All medium properties such as the constitutive parameters and 

the geometrical configuration are all absorbed into the reflection 

TM 7E* 

coefficients R and R , which is readily computed by Eqs. 
(2.31) and (2.32). Clearly, when all regions possess the same 
constitutive parameters, namely when there is no stratification. 


r™ = r te =o. 

2) Since the anisotropy in permittivity appears only in e (+) 
and the anisotropy in permeability appears only in p(+) , it is 

*■ ' *” ’ mjj 

seen from (2. 31)- (2. 34) that R does not depend on the magnetic 

TE 

anisotropy and R does not depend on the electric anisotropy. 

TM TE 

Both R and R are seen to be even functions of k . 

— . P 

3) It is obvious from (2 . 35) - (2 . 38) that a vertical electric dipole 


excites TM waves only and a vertical magnetic dipole excites 
TE Oaves only. Whereas both horizontal electric and magnetic 
dipoles excite both TM and TE waves . In the case when 
permeability is isotropic, the TM waves are extraordinary 
waves and the TE waves are ordinary waves. A turnstile 
antenna which consists of two dipoles perpendicular to each 
other and driven 90 degrees out of phase also excites both 
TM and TE waves . 
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4) The above formulation can be compared with the potential 
approach for the various cases that are existing. In the case 
of no stratified medium, the results are checked by using the 
identities 


(1) ’ + 1_ H (l) _ H (l) 

L k d 1 o 


(2.39) 


and 


i / dk e 2 (k p) = ex P (ifc/ P 2 + z 2 )/% 2 + z 2 

(2.40) 

2 2 2 

where r = p + z in spherical coordinates. We define the 
number of layers of a stratified medium equal to the number of 
boundaries. The medium below the nth boundary is called the 
tth 1 ay e r . 

5) The one layer case (or half space) has been extensively 
studied. We obtain from (2.31) and (2.32) 


R™ 


TE 

R 


e(+)o 

^<->0 


exp (i2k^ e ^d ) 
z u 

exp (i2k* m) d Q ) 


(2.41a) 


(2.41b) 


In the case of a perfectly conducting half space, 

TM TE 

e = c- — > co . Eq. (2.41) gives R =1 and R - -1. 
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6) Observe that all contributions due to all layers below the 
first are lumped into the reflection coefficient R n _^ * such 
that 



- exp(i2k z (e) d o ) 


(+>: 


e (-); 


(1 


e( + )o/e(-)o 2 e<->o/«<+>o . 

e(+)J/e(->J + exp(.i2k^ ) d 0 )R™ 1 


(2.42) 


# TE 

n 


exp (i2k, 


(m) 


d o> 


(l(+) 


li(-) 


(1 


»(+)><-)* - u<->>(+)* 

li(+)>(-)^ + exp(-i2k^'d 0 ) 


(2.43) 

The definitions for R™ x and R™ x follow directly from (2.31) 
and (2.32). 


In order to make effective use of the results obtained for 
the one layer case, we follow Wait (1953) to define a stratification 
factor for conductivity dominated media such as the earth. For 
a vertical electric dipole on the ground, the electromagnetic 
field quantities can be written in terms of one single integral 


00 



XfKo 1 ’ <V )e 


ik z 
z 


k z + Q e k lz 



(2.44) 


The field components are 

H . (- I* ) !£« 

H <(, . ' Wtt J 3p 


It 


3 2 J 


t ? _ / _ XJC ' } r. 

kp ' 5 n u e" 3p8z 


- <- . «=■ > £ i? * k2> J < 


(2.45) 

(2.46) 

(2.47) 
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The integral J follows from Eq. (2.35). The stratification 
factor. Q e is determined to be 

1 - R™! 

°e " . . -TM (2.48) 

1 + Vl 

Under the assumption that the stratified medium is conductivity 
dominated such that k{®* = k. t Q e is then seen to be a constant. 

1 Z 1 v 

k /k „ , and J becomes 

pz qz pf q e 


J 


e 


ikr -j -p, 

2 - — [1 + (irp^ 1 '^ e erfcfiv^)] 


(2.49) 


Eq. (2.49) is obtained from (2.44) for far field and observation point 
near the surface such that kp » 1 and 7/p << 1. In (2.49), 

r 

k p k. z 2 

P, * i — y- , P, = p. ( 1 + r (2.50) 

1 2k^ 1 o 

and erfc is the complementary error function. Replacing k^ 
by Qgk^ / we obtain an analytic solution for a vertical electric 
dipole on top of an N layer stratified media. 

For the case of a vertical magnetic dipole on a stratified 
medium, all electromagnetic field components can also be 
written in terms of one single integral . The stratification 

factor can similarly be determined from the reflection 
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coefficient 


r te 


Under the assumption of high conductivity. 



1 - R 


TE 


h -1 


1 + R' 


,TE 


h-1 


(2.51) 


m 


V O 11 ^ 0 P) elkzZ 

k z + k l 


(2.52) 


Eqs. (2.51) and (2,52) are dual to (2.48) and (2.49). The stra- 
tification factors for a two layer ground has been calculated 

by Wait (1953). In Figs. 2.1 and 2.2 the amplitude and the phase 

C> are plotted versas /wyo. d. for the two layer case. C> 
m lx m 

is seen to drop to near unity as /wycf^ d^ becomes larger than 

1. Thus the effect of the second layer diminishes as d^ > (wyo^) 

In Figs. 3-6, the amplitudes and phases of both and 0^ 

are depicted. The depth of the first layer is assumed to be 
—1 /2 1 /2 

d^ = 2 ( 0 ) 110 ^) . The ratio ' is set to 0.5, and 

1/2 

( 03 / 02 ) is used as a parameter. As d 2”^l becomes larger 

than 1.2 ( 0 ^ 03 ) , both |Q e j and |Q m | dro P to within 

about 6 percent of its final value. Beyond that we can 
safely omit contributions from the last layer. 

7) The integrals as presented in (2. 35) - (2 . 38) can be directly 
programmed with a computer or analyzed analytically . Asymptotic 
evaluations of the integrals and numerical results for the 
electromagnetic field quantities under various circumstances 
constitute topics of the subsequent papers. 


- 1/2 
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Chapter 3 . 

DIPOLE ANTENNAS ON A HALF SPACE MEDIUM 


3-1 FIELDS IN THE LOWER MEDIUM 


For a horizontal electric dipole above a lossy dielectric 
half space, the electromagnetic field components in the lower 
medium (i.e. the field components of the transmitted wave) are 
listed as follows: f 


— H.^ (k p) cos<J> 
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g** = 


/ dk 


I£o)y , 


8tt 


-ik 


(m) 

tz 


v t /v + k tz >/k z 


h] 1 ’ ' (k p) sin* 

r ^ r* 


(3.1) 


= / dk 


ii 

8tt 


V * k tz >/k z 


~ ik tz >z 

e sin<f> 


, (m) 

— ' (k^f ) sin 


(m) 

(k p) cos<(i 


~N 


k z p 


(3.2) 


k 2 

— £— (k p) sin<J) 

b k 1 p 
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E™ = / dk 


It 


p 8irwe 


ft + ftz 


_i y <e) 2 

2 ^ lk t2 2 

—XeT e 


- k kj® } H 1 (1) ' (k p) cos* } 
p tz 1 p I 


c (e) 

h| 1) (k p) sin* 


k. 


- - k 2 H* 1 * (k p) cos* 
a p 1 p 


0.3 


y 


r 




= / dk 


lte 1 

8-rre 




TleT 

L tz 


H l 1 ^ ^ k p p ^ sin( 5 > 


k H * 15 (k p) cos* 
pi p 


( 3.4 




A typical integral to be solved takes the following form 

I ft(kf)e il>nt H,‘ l> (kpf) diy } e< O 

The first order Hankel function is the negative of the derivative 
of the zeroth order Hankel function. In the above integral 




fe„ 

/«, + k„ 



Two points are to be noted: 

1. The path of integration along the real axis of ^ 
is actually undefined for strictly real k because 
the branch point k p ~ k then lies on the path. 
There are two remedies . One is to assume a vanish- 
ingly small imaginary part for k. The other is 

to choose the Sommer feld path of integration. 

2. Taking remedy into account, we notice that there are 

two double-valued functions in the integrand viz . 

*k and k. . They must be well-defined along the 
z iz 

path of integration. They are defined as follows: 
Along the path, /?c t Xat >0 j 

Ke k t >o * Ln >0 . 


- To facilitate solution of the integral 

the conformal transformation * k, si«/3 is made , so 
that the integral becomes 


I * k, J c sin/3 /I (ji) e ' ^ (k,/> si* i/s) 




where 

fe. cos/3 
* k x co‘i^k i 

C is the original path in the J3 
plane. If k^ is strictly real , C first runs from 
-J * i oo to -f- and then along the real /$ axis 
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to J and then from ^ to ^ - i oo . When is 

slightly conductive, C is slightly shifted with 

its position as shown in figure 3,1. 

Now, in this plane, only 

I -/IF k, 2 $i**/S is a double-valued function. The 

branch cut Im k_ = 0 is drawn. On the top sheet, 

z 

Im 0 , and on the lower sheet Im k 2 £ 0 

The path of integration C lies on the top sheet. The 

path C is then detoured to a path such that along the 

new path, the argument of the Hank el function is always 

big so that the Hankel function can be expanded asymtot- 

ically. 





* e 


L (k^sU/i - £0 



. m - i4 ( * cos/3 + 

S<a /i A (p) e 


( I 8 i ^ 


0 + 8ikfi\*p) dp 


The exponential arising from the 
asymtotic expansion of the Hankel function, combines 
with the exponential in front to give Lk, fccos (jZ~ 
where oL * • Therefore, there is a saddle 

point at . The branch point due to k z is at 

q* a qJ + i qJ 1 . Because k 1 has a small imaginary 
part, therefore is complex. For slightly . 

conductive material, B 0 * is approximately the critical 
angle. Q* , on the other hand, is very small but negative. 
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The location of the branch point and branch cut are all 
marked in figure 3.1. 

The problem, that is next faced, is to 
detour the original path C to the saddle point steepest 
descent path. There are two different detours for two 
different regions of interest. For smaller than “*> , 

the detoured path starts at - ^ + i , goes to -'J * + £ eo 

, and then goes down the path of steepest descent 
passing through the saddle point crossing the branch cut 
twice and crossing the real axis at the saddle point in 
between, and then goes down to ^ +* - too and 

from then back to the end of the original contour at 
^ -£eo .On crossing the branch cut once, the 

contour gets onto the lower sheet. On crossing the 
branch cut again it crosses back onto the upper sheet. 

The dotted portion of the path lies on the lower sheet. 

The integrand is analytic within the domain 
of the enclosed contour find along the path, and by 
Cauchy's integral theorem, the line integral of the 
original contour is equal to the line integral of the 
new contour. Notice that although the domain enclosed 
by the two contours lies on different sheets, the function 
is still analytic along the contour and in the domain . 

On jumping to the lower sheet, the function retains its 
continuity . On the other hand, if it jumps acoss the 
branch cut and stays on the same sheet, the function is 
discontinuous at the jump. Also, in this case, since 



3 


the new contour crosses the branch cut twice, going 
down and coining up again to the upper sheet, it can meet 
the other end of the original contour and forms a closed 
loop . 

i 

Therefore, for oC less than B 9 , 



where 

flw 


fc, COS * 

h,C6$ ot +7 


In the above expression, the double-valued function 

** Jk l - k* 5\r? << assumes a positive real part and 

a negative imaginary part. This is because the saddle 

point, as can be clearly seen from figure 2a, lies on 

the dotted portion of the curve and thus lies on the lower 

sheet. And , on the lower sheet. Ini k <0 . 

z 

Physically, the saddle point corresponds 
to the angle of observation in the lower medium. 

. _i 

For oC bigger than t? 0 , i.e. the 

observation angle greater than the critical angle, the 
branch point contribution has to be included and 
the original path C is detoured in the following way. 
(Figure 3.1b). It starts at * | t iw and goes to 


O + { DO 


goes down the left side of the branch cut 
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on the upper Pieman sheet, goes around the branch point, 
back to O + too on the right side of the branch cut 
on the upper sheet. It then crosses the branch cut onto 
the lower Riexnan surface, goes to — ^ + * +• t and 
from there, follows the steepest descent path, crosses 
the branch cut again back to* the upper Rieman surface, 
crosses the real axis at saddle point ©C. , goes to ~ i 

and finally connects back to the end of the original 
contour at ” L oo . Therefore, summing up , there 

are two contributions to I : the saddle point path 
contribution 1^ , and the branch cut contribution I 2 . 

The saddle point lies on the upper sheet. Therefore 

! I,= f - [ * («> + ikz .+ cat* fl V“0)J 

where , 

k, Co JX 
Coset +7* 4 

and k assumes a positive imaginary part and a corr- 
z 

esponding negative imaginary part in the above expression 
since on the upper sheet, Im k )0 . 

At this point, it is convenient to adopt the 
following convention. We denote ( the 

square root sign ) as always having a positive real part. 
Therefore, going back to our analysis, k can be either 
+' Jk X -k l z siS« / or - depending on the sheet 

the path is lying on. VTith this convention, the solution 
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is rewritten as follows: 

i 

for cC less than , 


I * “T" f Af*0 + 2lfc t K (A'W + Cat* A'CO)] 


where, 


k t c os<* 


for ©t greater than 0* , 


iJ*,R f 

I, = T *r ffl (*> + 3 Tm + cat* 


» (*; * ■ 

fe.cos ec - y^ -Vsin V 


To evaluate the branch cut contribution 

I 2 # much care has to be exercised. The branch cut is 

defined by Im k = 0 . On the left side of the branch cut 

z 

Re k z > 0 , while on the right side of the branch cut 

Re k 2 ^0. Denoting the positive side by B + and the 
negative side by B_, 


nz 

* J ir/» 


-ItP f jk tk, Rcos(ft-ot) 

e Is t» (l fi(fi) e <ip 

aow« 8+ 

Up 8. 




ft Cos C/3 “ 


ife, ft tos C/3 -*0 ^ j 



( Up means from branch point to ©o * 5 

The above integral can be combined in two 
forms, either by chan going B_ into B + so that both integrals 
are written as B + , or alternatively, by converting B + 
into B to have both integrals in terms of B_ • Thus, 

there are two forms of ^2 : 


r. -yf e' if j s „ 

up 8* 


d/s 


or 




€ 


■if f i 
I 5* A 
1 


-2k, tk. i*,«e «(/>-«) 

E ■*/> 


Both give the same answer. 

The next step is to detour the path, either 
B + or B_ , to the steepest descent path passing through the 
branch point. ( Note: one must distinguish between the 
steepest descent path passing through the branch point and 
the steepest descent path passing through the saddle point.) 
The position of the steepest descent path BS through the 
branch point is clearly depicted in figure 2b. It lies to 
the left of the branch cut. Therefore B_ can be detoured 
to BS by simple application of Cauchy's theorem. Therefore 



S “'/* fe ' Co V 3e d P 


On the other hand, if B + is to be detoured to BS , one needs 
to cross onto the lower sheet and subsequent calculation 
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needs to be performed on the lower sheet. This is unadvisable . 
If t it is quite plausible, one is not careful enough to 
follow the above step by step, one can easily just detour 
B + to BS and arrive at 






2k, 





This is the negative of the correct answer and is wrong. 

The contribution of BS mainly comes from 
around the branch point. By applying the method of 
steepest descent to the above integral, one obtains (Appendix 


I) 


r. ■ 


+ ikf 


[>' ** 0- Cot J. tan B t ) yz 


- A, 1 - k‘ 


Finally, in summary, 
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for cL less than 
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[am + jik$ 


(A*M + cotec 


where 

kj co s e c 

k,C0S< +Jk i -k l x SmV 


for oC greater & 0 



Im k t <0 



3-11 


I * “f* 6 ft fd + a/fe.R (A*C«0 * Cot ec /T(*))J 

/ . 

* JV e' iM * U/ * _ | 

• ~^ r <x ( i ~ cat cl to* 6J'* j IM kg 

where, 

fe,C6$ol 

Au > ’ k,a 

Several physical observations can be made. 

For oC greater than d 0 * , there are two waves, repre- 

sented respectively by the first and term in I. The first 
term corresponds to the spherical direct wave, and the second 
term corresponds to the flank wave . From the equation, one 
also notices that if k^ is strictly real, both terms blow 
up at e(. - 9 C . Mathematically, this is to be expected 
because the saddle point collides with the branch point. 
Actually when the saddle point approaches the branch point 
and gets near to it ( this also applies to the case of 
complex k^ where oC is never exactly equal to B© .) 

, the saddle point method breaks down. The area cC ^ B 0 
, is the area in which the two* waves interact and they have 

to be considered together and not separately. A detailed 

3 

discussion is given in the book by Brekhovskih 
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by most authors. The advantage of using this branch cut 
is that it can be detoured to the steepest descent path 
passing through the branch point. In the following# it 
is shown that the same answer \^ill be arrived at even if 
a different branch cut is chosen. The different branch 
cut is Re k„ = 0 . Solution with a different branch cut 
yields the same answer means that the choice of different 

/ i 

signs for k_ in the two regions cC < B a and oC > 6 0 
is absolute and is not arbitrary. 

With the branch cut Re k = 0 , and after 

Z 

the = k^ sin ft transformation is made# the same integ- 
ral is obtained but with the picture in the complex 
plane changed as in figure 3.2. BS is the steepest 
descent path passing through the branch point. On the 
top sheet Re k > 0 # and on the lower sheet Re k <0 . 

On the top sheet, where the original path lies, there is 
the positive side of the branch cut B + where Im ^ 0 # 
and the negative side of the branch cut where 

In k z ^ 0 . The most important observation is that, in 
this case, the branch cut , no matter whether the B + 
side or the B_^ side, cannot be detoured to the steepest 
*■ descent path BS through the branch point. 

Analysing the integral as the case with the 

branch cut In k = 0# for ed < 6L , the saddle point 

z 

lies on the upper sheet and Re k )0 . 

z 

I “ T — ( fl (*) Ji fe,R + cot oc A'(cO)] 
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Notice that the saddle point contribution in this case 

i 

is different froin the corresponding 1^ in the case 

In k =* 0 because in that case Im k ^0 and in this 
z z 

case Re k > 0 . 

z 

Analysing I? the branch cut contribution. 


i, = 


rar f . t . , (p-<a , , 

m irp & J S‘n pi fl (p) e ( I mkfpsinfi) <(p 


down 3. 
uf 0 + 


rzhT ~ l * f . Cte.R tosffi-oc) _J 

6 J sin P k t co^p e C I + 


r 


uif> 8. 


2 k, 




Df/i) 


( Again, up means from branch point to do . ) 
where , 


2 fe* k, 


* fe.W/3 - k* + fc, z s.V/i 


tfe k t > 0 
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whe re , 


k, CCS u, 

A ~ k, cos oc +■ STO 


i k £ >0 


Therefore, the answer checks with the case with branch cut 

In k = 0 . 

z 

In the case ^ > bJ , there are , as in the 

case with branch cut In k = 0 , two contributions ? one due 

z 

to the steepest descent path passing through the saddle 

point and the other due to the branch cut. The branch 

cut contribution, however, as opposed to the case of branch 

cut In k s o , is of first order effect, 
z 


then , 


where 


Set 

1 = *1 + *2 

1^ = contribution due to steepest descent 
path through saddle point. 

Ij = branch cut contribution 



[flM + jTk^ ( + Cot * 


where. 


A( -- 


k, COS « c 


t0 J * Vh v - kJTUTTZ ■> 


R. 


K >o. 
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The path B cannot be detoured to BS . To evaluate the 

i 

contour integral along B_ , the path is detoured to 

a path consisting of two parts: the first part extend- 

ing from the branch point to the saddle point, and the 
second part extending from the saddle point to 
down the steepest descent path. By Cauchy's theorem 


l "1 + l- 


Therefore, 



-1 


+ J 


X r £ ® 2ik t R, ( & (*£■) +• Cot o(. 

~t«Z -f-t/y* 


~h 


_ e 

a 


f* z «».* ( l - cob oc ta* )*/£ 


where 


E>M 


£ /fe* - fe, * S»»* fe , C05 *£ 

— k* + A,* Sin* oc 
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Combining these expressions, one has 
* 

i - i, h 


i. A 

i 


ik ' R r I 

H ( dr (°0 +• 2 ik,R ( £ + 6 * (*) Cot « 


„ . -to* +- ikf> 

i 2 w. e 

|> a a* (I - cof eC ton 9 0 ) $/2 



where Q(^) ~ /) (tl) +• D («0 


A ( oC ) of the case of branch cut Im k = 0 . 

z 


Therefore/ in summary, the answers are identical 

whether the branch cut Re k = 0 or Im k = 0 is 

z z 

chosen'. 


For TM waves of nonmagnetic material, 

te- 

on the top sheet, there is a Sommer feld pole at sih ~ 
in the denominator. Its effect can be neglected by the 
following reasoning. VIhen U. , K, &<» t ^ lies on the lower 
sheet { branch cut Ira k z = 0 is chosen. ) , while lies 

on the upper sheet. Thus the pole is not enclosed between 
tKe original contour and the detoured contour. The pole 
does not contribute anything. VJhen oC > 0 O , the 
analysis is slightly more complicated. When the original 
path C is detoured to the saddle point path p , two 
contributions need to be considered : one due to the branch 
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cut B, the other due to the residue of the pole P 1 be- 
cause in this case, the pole does lie between the original 
contour and the detoured contour. But, in evaluating B, it 
is detoured to the steepest descent path BS, and on this 
detour, P^ is again enclosed, and this time, its contri- 
bution is exactly the negative of the residue and thus can- 
cels out exactly the previous contribution. Summing up, P^ 

does not contribute both for a < 0 Q ' and a > 0 Q 1 . The 
only contributions come from the saddle point and .the branch, 
both of which have been discussed already. 

3.2 RADIATION PATTERNS 

In this section, we analyse the radiation pattern 
of the horizontal dipole for the half space case. The radi- 
ation pattern gives information of the angular distribution 
of the radiated energy. This information is of crucial im- 
portance if the interference patterns for the two layer case 
are to be explained in terms of ray optics. 

The radiation fields are of 1/r dependence. The 
1/r term, in this case, corresponds to the first order term 
in the saddle point contribution. 

For the half space case, with the dipole on the 
surface, the reflection coefficients are 

hf _ n t k t - ufe, z . o™ - C, kt ~ Sfe, z 
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The electromagnetic fields , to first 

order are the following. 

For the upper medium: 


TE 


_ ~ -2 (, k sin ? X 0| (*) Sin 0 

4 w r 


H f - 


T£ 


= M * * 


H K s “4nT "V - ^ C0S 6 X oi ^ S,>1 ^ 




- 



4-tt t 


- IL _€ 

ss 

4ni 


1/ 


” 4niw£ 


If . 

— 

4-rfc>£ 


IL € 

rr 

4iri 


i/*r 


; 4 »- 


TM 


r* _ _ 


TM 


— K 1 cos 1 B Y lo (9) cos 0 

i V 

— ik* $•* & cos 6 Y lo (&) cos 0 

ikh 


For the lower medium : 


H 


T£ 


_ IL e ik,r • / ■ V , N . . 

4tr r L Sur\oC A J0 £«0 sm 0 


T£ It 






TFt 7 X l0 SIM 0 


. . TE It p , \/ , \ 

Hp * TFT “T“ fe. cosod l 0 ^ Sin 0 


J. ™ _ - I*. . « ,1 * v , , 

Cp 4-ITituSi r «•» CoS od | #( ( tO co$ 0 


TM II o ii *< r 


p * s c * / 2 , i* 

k* ® k-TTMB, f L *« Sl>t^COS«L } O/ C°c)cos 0 
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H 


7M 


11 

4jri 


ik,r 


k, CO 5 oC Y ol C<<) C*s 0 


where r 


fa; 




,v lJ 

jft fe i2 

+A 

k j* 

Sy 

35 -^i— U 

- Si 

+ 6c 

hi 

k J* 

Xy 

= / + 1 

*0 

* - 




2 £; 


= 1 + 

$cj. 

- J 


0 is the observation angle in the upper medium and 
o L is the observation angle in the lower medium. 

The Poynting vector is s = Red (EXH* 
Denoting the radial component of the Poynting vector by 
S r , we obtain , for the upper medium 

s r •tigffektiw «*’«• +• fc c* l e r>«v j 


for the lower medium. 



"2 (W (T^ 






3-20 


The maximum power radiated into the 
lower medium is at the critical angle for broadside 
( <P ® *10° ) and at the angle 

for endfire ( 0 * 0 * ) . These values can be obtained 

by setting the derivative of S r with respect to the 
observation angle equal to zero. Another interesting 
observation is that for the endfire radiation pattern, S r 
is precisely null at the critical angle. A quick way to 
explain this is that the free space radiation power is 
zero along the endfire surface, so that there is no 
" incident field " into the dielectric at 0 * SO* 
when the dipole antenna is put on top of the dielectric 
medium. Another observation is that along the surface, 
fields are of second order effect, and thus the radiation 
power is zero along the surface. If k^ contains an 
imaginary part ( for conductive medium ) , however small, 
the radiation is hull in the lower half plane, because, 
by definition , the radiation pattern is measured at infinity. 
The Im gives rise to an exponential attenuation with 

distance and the power vanishes at infinity. 

3.3 FIELDS IN THE UPPER SPACE 

In actual experiment, the receiver is 
placed above ground, so that, we are actually interested 
in the region z ^ 0 for the half space problem. The 
transmitted wave for z ( 0 is interested not in the 



3-21 


half space problem, but, as has been mentioned, its 

solution sheds light on the two layer problem which v;ill 
> 

be discussed in the next section . 

For z ^ 0 , a typical integral to be 


solved is 




(k rP ) 





Notice that this integral is different from the trans- 
mitted wave integral in that k, is replaced by k 

lZ 2 


Let 


-TThc 


In this case, the conformal transformation is 


kp = k 

The saddle point is at 9 


where 


9 


is the obser- 


vation angle in the upper medium, 0 - f- ah ~ 

The branch point is at 9, ~ $ in' -g- = Q ' - i. 

The branch cut chosen is Im ^ = 0 . The equation of 

the steepest descent path p passing through the saddle 
point is /?* feR cos (ft- 6) * feR . Let’s take ^ to 
be strictly real for the tine being, so that q/ = X 

cosh 9, ~ . The* steepest descent path passes 

through the branch point when 

Re kR cos (&,-&) - kR 

cos (9/ - 9) cosk 9 t a = 1 
si* 9 = 
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and this turns out to be 0 o , rruite a coincidence. 
Therefore, for 0 < B c t only saddle point contribution 
is needed, while for 9 > B 0 , branch cut contribution 
needs to be included too. The branch cut contribution 
is to be carefully manipulated by the B + and B_ 
approach, and also by using the method of steepest descent. 
For the Tm wave, the Sommerfeld pole does not contribute 
by similar reasoning. 

In summary, extending also to slightly 
conductive tOO/ 

for 9 < & r , 

I = T ^ f B (9) + 2 ^* ( &'(9) c at & B"(B J)j 


for e >e* 


o ikR. 


r = 2 e 

1 i R 

[B (8) + 2ih.R ( 8 (8) cot & + B" ( 0}jj 

2 k, 

+ ih,p 

e 

1 

h> 

* 

p 

At 

(i -hi sec B 0 c ot$)¥2 


The double-valued function k^ z has its imaginary part 
bigger than 0 on the upper sheet. Since for both regions 
of interest 9 < 6 0 and 9 > 9 0 the saddle 

point occurs on the upper sheet, therefore Im k^ z 0 



in the above expressions. 

As can be seen from the expressions 
for / there are two waves , the first term re 

presenting the direct wave , and the second tern the 
inhomegeneous wave. The inhonegeneous wave decays 
exponentially away from the surface. Its magnitude 
is greatest on the surface where it interferes most 
significantly with the direct wave . To do detection 
of electromagnetic properties, we know that for two 
waves with different wave vectors k and k^ , the 
peaks and troughs are separated by a distance of 

* * HvT 

On Lhe Ollier hand, if we know , 


then k^ can be determined. 

On the surface. 


T - _ 2 k, 

Actually, the exact ans\-;er can be 
obtained by the following steps of algebra due to 
Van Der Pol. Use is made Qf the following identity 


. /■» 
jr -*/ • 


//) 

If- e 


The integral is 


1 = Lx^ H ° (h rr ) dk r 
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r CO 

~ V-k 1 k ti ) H 0 > ( kp p ) dkp 

C k r k ‘^ ( ^ dk r-[-£i^l 0 


L V'^o^uv * [~h f nr-J f , 0 


•'• 1 * £ [T ^ (^-p-) ~ T- S ^( ik '~t>] 


The various interference and radiation patterns due to different 
half-space medium compositions are presented in Appendix A. 
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FIGURE CAPTIONS 


3 j a. The complex |3 plane for the transmitted wave 
with the observation angle smaller than the 
critical angle . 

b. The complex ^3 plane for the transmitted wave 
with the observation angle greater than the 
critical angle . 

3 2 The complex |3 plane for the transmitted wave with 
a different branch cut. Re k 2 = 0, is the 

contour joining the branch point to the saddle point, 
Lj is that part of the saddle point steepest descent 
path extending from the saddle point to -i«* . 

V is that part of the saddle point steepest 
descent path extending from +i«o to the saddle 
point. 

3,3 The complex p planes waves in the upper medium in 
the half space case. The transformation is 

«= k sin p . The branch cut is Im k l2 * 0 . 









4-1 


Chapter 4 . 

GEOMETRICAL OPTICS APPROACH 


4.1 TWO LAYER MEDIA 

We first consider solutions to the reflection coefficients 
for a two layer medium: 




V Roi ~ ^Qi 


yR ot +-R u «x / »(i2^ li «t; 


-i 


and 


R 


tm | 




r t - ~ 1 

** ^ V^oi + $1Z &*P t ^ ^1* ^0 J 


1 + 


0 tE 1± ft? e/p(t2k 1£ cO v 

^ ~ i f ^ OI ^ 


e*p(i 


i-/? 


TM 


l-5„eyp(t?fe lg J-) y 
i) ' io 


In the geometrical optics approximation, 
the interference pattern is determined by three waves: 

1. the direct wave from the antenna 

2 , the inhomogeneous wave excited by the antenna on the 
surface, and 

# the wave reflected in the subsurface following ray 
optics paths (Figure 4.1). 


3 
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TE 

The denominators in 1 + R and 
TM 

1 + R are expanded in power series : 

1 f /? TE « X M [ u- X„C m<i) ] 

1^-R™ = Y l0 ( 1 ~ ^ T ol S 0 S 2 exp(i?k u m<t)] 

Ohe first term in the power series gives rise to the half 
space solution which consists of the first two kind of 
waves, both of which have been dealt with in the section 
on half space. Mathematically, the direct wave corresponds 
to the saddle point contribution, and the inhomogeneous wave 
correponds to the branch point contribution due to k^ z . 

The summation term gives rise to all reflections from 
the subsurface. At each observation point, the number of 
reflections from the bottom layer can be determined by 
exercising ray optics. Each optical path can be traced 
back to an image source . 

The integral for the "reflected 
terms of images" is of the following form 


h ** H> V^p)X OI X ia e* a,, '’ u< K" o 1 r2 


By rearranging terms in the integral 


r 5 = I, [ dk r-t { h D ( V X « K. C C e 
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If z is small, the important observation is that each 
of the integrals in the summation is of the form of 
the half space transmitted wave, ( that is why the half 
space transmitted wave is solved first.) , where 2md 
now plays the role that was played by z in the half 
space transmitted wave case . The ~ k, 5w»/3 trans 

f 9 rmation is made, the saddle point contribution is. 


I, f r 4^-0(k r ) ^ X l0 OC e ‘ V 


tn*l 


S^oC n * 


R-m ~ tZndif 

" +an '‘£r 

The saddle point is at - tan 

The above results hold for k^d » 1 and d»z. When 
exceeds the critical angle, then branch cut 
contribution needs to be included. The branch cut is 
evaluated by deforming to a steepest descent path through 
the branch point. In geometrical optics approximation, 
integrals along the BS's , the steepest descent paths 
passing through the branch points , are usually neglected. 
In deforming the branch cuts to the steepest descent path 
through the branch point and before deciding that the final 
contribution due to the latter is negligible, it must be 
shown that the deformation is valid. If the deformation is 
invalid, two alternative approaches can be taken. One 
alternative is to choose another path and the other 
alternative is to choose a different branch cut. 
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In the complex (J plane, there are two 
double-valued functions, k z and k 2z , and their values 
are to be defined. In analogy to the transmitted wave 
case, therefore, for 

-s. 

1m. K 0 5 

for > Qj 

In 

By the same argument as presented in the 
transmitted wave case, the contribution of the Sommerfeld pole 

fe. 

P^, , ( Sifl * ^=f==|= ) , is well taken care of by the 
steepest descent path BS • The net contribution is 
the saddle point path p and the branch cut steepest 
descent path BS . BS is to be neglected in geometrical 
optics approximation . 

For the double-valued function k 2z , 
it is to be discussed under two different cases: one 
case with a good subsurface conductor, the other case 
when the subsurface is a good dielectric, 

Whpn the subsurface is a aood conductor, 

. -1 ^2 

the branch point 0 2 = sin £— and various singularities are 
shown in figure 4.2. For this case, the branch cut Im k 2z = 0 
is chosen because this branch cut can be detoured to BS 2 , the 
steepest descent path passing through the branch 0 2 . BS 2 
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is of second order effect so that it is neglected in 

geometrical optics approximation. With respect to ^22 # 

there are two sheets. On the ‘first sheet, on which the 
* ^ 

original contour C lies, Im k 2z ^0 * on t ^ ie bottom s ^ ee ^ 
Im k- < 0 . Since the saddle point occurs on the first 

4* Z 

sheet, therefore the value to be taken for k 2z * * n 

above expression for saddle point contribution, is such that 

its imaginary part is bigger than 0. The Sommerfeld pole 

P 2 , ( Sin m 5 ' aS Can be S6en fr ° m itS 

position in figure 6, is never included between the 
original contour and the detoured contour. It, however, 
may affect the saddle point cd m if is getting close to 
*^p , In our regions of interest of big d and 

moderately large p , this is never the case. Thus, the 
effect of P 2 can be safely disregarded. Therefore, 
in suirmary, the values of double -valued functions k z 
and k„ are to be chosen in the following manner. 

for oCm < 6 0 > Im kg < 0 , Im k 2i >0 

for > 0„ , Im kg. >0 i Im kjg y 0 

When the subsurface is a good dielectric medium, the 
various singularites and contours are shown in figure 4.3. 

The following observations are made: 

1. This case contrasts with the previous case in that 

- -* h 

0, * $tn -r^* is at a drastically different 
2 *1 

location . 

2. The branch cut Im z ~ 0 cannot be detoured to the 


1 


V*f 
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f 

steepest descent path BSj because the detour involves 
passing through regions where the exponential blows up. 
It is to be emphasized that, in geometrical optics 
approximation, it is the BS ? , steepest descent path 
passing through branch point, that is negligible, and 
not any arbitrary branch cut one defines. The pre- 
scription, in this case, is to define a different 
branch cut such that it can be detoured to BS 2 . The 
branch cut chosen is Ee k 2z = 0 which can fulfill 
the function.. 

3. With respect to k 2z , again, there are two sheets. 

On the first sheet Re k 2z >0 » while on the second 

sheet. Re k ^ 0 . The branch cut Re k, = 0 

can be detoured to the steepest descent path BS 2 
which is neglected in geometrical optics approximation. 
Since the saddle point alv;ays occur on the upper sheet 
with respect to this branch cut, therefore Re k, ^ 0 
is the value to be taken in our expression. In summary, 
the values of the double-valued functions k and k~„ 
are to be chosen in the following manner: 

for , Ir* < 0 J Re k Jt > 0 ; 

^ or ^ 0* , ^ )0 , Re ^0 
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The two layer case is different from the 
half space case in one other important aspect. In the 
half space case, the receiver is close to the surface, so 

JT 

that the saddle point is always close to 2 .On the 
other hand, for the two layer case , the saddle point for 
the "reflected terms of images" have their saddle points 
moving along the real ^3 axis as ^ and 
2md are changed. When the saddle point approaches the 
branch point 9^ , the saddle point method crumbles. Using 

Watson's Lemma, the criterion for the validity of the saddle 
point method is found to be 

j fc, (sib &o -sin oQ*| ^ 

The above condition can be achieved by 

either making k^ more conductive so that has a bigger 

imaginary part or by making R m bigger. But both of 

these procedures will lead to a decay of the reflections and 

thus decrease the interference effects; they are largely 

% 

undesirable. 

4.2 SOLUTIONS FOR THE FIELD COMPONENTS 


By making full use of the mathematics 
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developed in Appendix I and previous sections, the field 
components have been derived „ In the vicnity of the 
surface, i. e. z«<£, p , the six electromagnetic field 
components are listed in the following: 


H™ * fsm e (ik - X n (W + b c°$ex' ol (e) 


2r X„ (6) J sin 0 


IL Jt 2*,* £»+** . w 

^ f* 0 "* C 1 1 i Sin ooi $) ^ ' 










v* 
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E> t£ = J *,T V“ f X n (6) (<*t s X ol (e) +Xo t ( & J)J sinfi 


. Uum 2k, £^jV „ . 

9-jr /> ! »‘ (IH st c B„ccte) V2 ^ S,n ^ 


+ K CCe^l 




#n=l 
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r“ l ?(&> + J7*; (cotet>'(6) + f* 0 (oy)] ±™<p 

XL i£. e - ai ^ ik '/ > . 

~ fl />v ®; 3/1 sw '^ 


- /* ft y y c> m * o K **«* 7 : ^ 
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t -TM 

F* ' <Hrw£ 


[ (ife “TO Slii 6 0 (0) + " > “ o' f 0J 


.+ ItF s ‘ tt ^ QVs)] cos^l 
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71+ is«c e. cm ejv* 
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E™ = + adTr Uot 9 a; 16) +Q l "u»)Jcoi? 


It 

■J* 4ttw£ 




^-a* ♦* i^xp 
TIT i sec % cot e) v* 


cos 0 





fej 6W «t*. 


,ik,R* 


COS 0 


ty,™ - "^f -^T— f a (9) t 5 s; ( cot 9 Q '(«> + Q* (9) )j «s i» 


ii zi k * ± 

4*1T y*o- { 1 +■ * sec 0» cot OJ 3 ^ £i 



CoS 0 


whe^t P(0) - k CoS 0 X 01 (0) 
Q(9) = cos e Y u (S) 

Q,(9) * k* *** » Y»(») 



4-11 


Seyeral points are to be noted in the above fomula: 

i 

1. The first term is the direct wave. The second term 
correponds to the inhomogeneous wave . The third 
term correponds to the reflected terms of the images. 

2. In the first term, which is the saddle point contribut- 
ion, the first order term correponds to the radiation 
field. 

3. At z , i.e. the receiver being close to the 

surface, the terms 5 ^(B) and Q (B) are much smaller 
than their derivatives , X<, A ( 0 ) t X 01 (&) * Q^B) / and &"(&) . 

This is because of the cos $ dependence. 


4 . 3 NUMERICAL RESULTS 


The interference patterns for the six electromagnetic field 
components for different configuration and media are illustrated 
in Appendix B. In the diagrams, the magnetic field components 
have been normalized by 


n - 11 

H ~ 4irV 


u 


while the electric field components have been normalized by 



i F 

4- V A 


The interference patterns may be best 
understood by separating the contributions of the various waves. 
Thus we depict in Figures 4.4 - 4.5, two typical field 
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components H and E separated into the half space 

solutions and the waves reflected from the first boundary . 

We note that ll 2 together with Up and E^ are the field 

components measured broadside, i i .e . 0 ~ 'j? • while 

, E p and are those measured endfire, i .e . 

^ = 0 . The maximum of the first reflected waves for 

H and E occur at different locations, 
z z 

To explain this difference, we consider 

the radiation pattern due to the dipole lying on the 

lunar surface. Under the geometrical optics approximation, 

the maximum of the wave reflected from the subsurface 

cor reponds to the angle of the Poynting power that is 

radiated into the first layer. For a HED at <!> = v/2, the 

maximum power occurs at the critical angle sin ^(k/k^) 

and k^ is assumed to be real. But for ^ =0, we find 

-1 / 2 2 ~ 

that the maximum -occurs at an angle sin J 2n /(I + n ) , 
where n = k/k^ 

The following observations are made: 

1. In the direction of maximum power, not all field compon- 
ents have to be maximum, e.g H z and Up have max- 
imum at different locations. 

2. A change of permeability from to^u. merely 

changes the magnitude of k^ . If jx, becomes 
complex, then the positions of the various singular- 
ities change, and the integrals must be re-evaluated, 
on the new complex plane . 



3. A small loss tangent change has small effect on the 


positions of peaks and troughs but changes the ampli- 
tudes of the field components drastically. 

The results of this thesis are now 

3 

compared with previous findings . Brekhovskih , in his 
solution made the kind of mistake due to the wrong 
approach to B + and B_ as discussed in the section 
on half space solution. He got the negative of the 
correct answer. Thus, if one plugs in z = 0 in 
equations 19.36 and 21.14, one obtains 

-L Jz. p • V . 2 

e + e 


and thus, fails to reproduce Van Der Pol’s exact answer. 

The validity of the calculations of this 

thesis can be verified by carrying out exact solutions 
of the integrals. Tsang (1972) obtained such solutions by 
numerical integration with a computer. The computation 
is extremely time consuming and a single comparison is 
shown. In Figure 4.6 , the exact solution may be compared 
with the .asymtotic solution under the geometrical optics 
approximation . It can be seen that the geometrical 
optics approximation is not good near peaks in the 
interference curves. Mathematically, the peaks are 
those regions where the saddle point cC m approaches the 
branch point 0 O . The saddle point method crumbles 

One more point requires discussion. 

In the geometrical optics approximation, BS and BS 2 , 



4-14 


the branch point steepest descent paths contributions, 
have been neglected. These contributions correspond 
physically to the reflections of the lateral waves 
from the bottom layer. Thus, it is implicitly assumed 
that the energy arriving at the upper surface of 
layer 1 due to the reflected lateral waves is small. 

For a better approximation, these effects should be 
included . 

The calculated results using the geometrical optics ap- 
proach are presented in Appendix B. Both a vertical magnetic 
dipole and a horizontal electric dipole are considered. 

Rossiter et al (1972), has performed the model tank ex- 
periment with a dipole over a tank of dielectric over a perfect 
conductor. The theoretical results are compared with experi- 
mental findings in Figure 4.7. 

Experiments have been done on glaciers 
to check both theory and experiment. The Athabasca 
glacier, located about 75 miles south of Jasper, 

Alberta, has been thoroughly studied by others with 
such methods as seismology, gravity, resistivity 
measurements and drill holes. In Figure 4.8 is a 
single comparison. The thickness of the ice determined 
by our theoretical curves is 180 meters which agrees 
very well with the thickness determined independently 
from the drill hole. 
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4.4 STRATIFIED MEDIUM WITH MORE THAN TWO LAYERS 

Cases of more than two layers can be 

dealt with in similar ways, using geometrical optics 

approximation. The steps are the same: the denominators 
TE TM 

of 1 + R and 1 - R are expanded and combined with 
other terms in the integrand. A suitable transformation 
is made e.g. c 

etc., depending on the predominant term in the exponential. 

By expanding the denominators in the 
n-layer reflection coefficient, one can get, in general, 
two terms: the first being exactly the same as the 
n - 1 layer solution, the second term corresponding 
to the reflections generated by an additional boundary. 

This, indeed, has been the case for the two layer problem 
where we get back the half space solution plus the "reflect- 
ed terms of images". The three layer problem is illustrated 
in the following: 
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.Thus, there are three factors, viz., 

v • x 

A 01 > ”, A io •'i© 9 c > 




Uz X„ X„ Ce Si ^' C ; f? 


w*a t*i P«» 


•x n? s 


(~tf {«»!?„ e 21 ^^} 


J>+t 


Th$e first two terms combined is just the two layer solution. 
To find the location of the saddle point in the last term, 
we observe that the exponential is 


i hpp + 2 i k ti m c/j +■ 2i kjj (-pit) ( d a - d x ) 

so that the two equations to be solved are 
kp ~ k, Stn 6 i ~ k 2 Stin 0 a 


2mJ l fan 0 X 4- 2 + (d A ~d J ) tdn B 2 ~ yO 


where is the saddle point in the - ^Si'yS trans- 

formation, B z is the saddle point in the ty - yQ 
. transformation. 

Note that physically, the above equations 
correspond to the ray paths given rise by the multiple 
reflections in geometrical optics when there are several 
boundaries of reflection. 
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For complex* k^ and k 2 , there are 
actually four equations and four unknowns as both the real 
and imaginary parts of the equations are to be taken • For 
slightly conductive k^ and k 2 , the real parts are solved 
first to get an approximate answer first. The imaginary 
pkrts are then taken into account and the final answers 
are obtained by perturbing around the approximate answers 


using Newton’s method of solving roots of equation. 



or the transformation 


bp = b, sw/3 


is made, depending 


on which layer is thicker. Care has also to be exercised 


in choosing the values of double -valued functions, the 
k z 's , in the complex plane as in the two layer case. 



FIGURE CAPTIONS 

4.1 The image diagram for the two layer case. 

4.2 The complex $ plane for the two layer case with a good 
subsurface conductor. 

4.3 The complex f3 plane for the two layer case with a sub- 
surface dielectric. 

4.4 The half space solution for H z and also the 
first few reflected terms, 

4.5 The half space solution for E z and also the first 
few reflected terms. Comparing figure 12 and 13 
show that the maximum for the reflected terms 

occurs at different location for the two components. 

4 . 6 Comparison of geometrical optics approximation 
with exact numerical integration. 

4.7 Comparison of model tank results with geometrical 
optics approximation. Dotted curve is G. 0. A. 
(Geometrical Optics Approximation) . Solid curve 
is experimental measurement. Vertical scale is 15 
db. per inch. 

4.8 Comparison of G, 0. A, and Athabasca Glacier data 
measurement. The stratified diagram gives the 
parameters of the G. 0. A. curve. Vertical scale is 
8 db. per inch. 
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Chapter 5 . 

MODE APPROACH 

In this chapter , we use the residue series method in combi- 
nation with either the modified saddle point method or the numeri- 
cal method to evaluate the integrals. Each term in the residue 

expansion can be attributed to a normal mode of the layered medium. 
The characteristics of mode solutions are then discussed. VJhen 
the non-modal part is evaluated by the modified saddle point 
method, the results are applicable to the far field regions only. 
VJith a combination of the residue approach and the numerical 
method, the results can be extended to all regions of interest. 

The advantages and the disadvantages of the various limits are 
also discussed. Theoretical results are illustrated for the > 
various field components plotted in the form of interference 
patterns . 


5.1 METHODS OF APPROACH 


Consider the integral 


H „ a k (l + R TE )e ik ^ "l 1 * 


(5.1) 


where H denotes the z-conponent of the magnetic field vector. 
VJe use cylindrical coordinate system with z-axis perpendicular 
to the surface of the stratified medium and P denotes the 
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transverse distance with origin at the transmitting antenna, the 
azimuthal angle <}> is measured with respect to the transmitting 
antenna axis. The propagation vector k in the free space re- 
gion has components k p and k z . I A denotes the antenna 
strength. (k p p) is the first order Hankel function of first 

kind. P. TE represents the reflection coefficient for TE wave 
components due to the stratified medium. 

TE 

For the two layer case, the reflection coefficient R 
takes the following form 

1 " (R oi >2 

r te „ 1 [! „ ] (5.2) 

R $f 1 + R Qi R i2 exp(i2k lz d) 

iii L> gip t 

where R^ and are the usual Fresnel reflection coefficients 

for TE waves at interfaces between the first and the second 
boundaries. 

In order to solve the integral in (5.1) , we make the trans- 
formation 


k p = ksinO (5.3) 

so that the two Riemann sheets due to the branch point k z is 

unfolded, l'7e then have four Riemann sheets due to 

Re k, and Re k„ . In the integral (5.1), we can write 
1 z 2 z 

H, (k p) = H. (k p)e~ lk P P and insert an exponential factor 

1 p 1 p 

exp(ik p p) which is in fact the factor for the I-lankel function 
at its asymptotic limits when k p p tends to infinity. Correspond- 
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ing to this factor, there is a saddle point at 0 = ir/2 . We then 

detour, on the complex 0-plane, the original path of integration 

to the steepest decent path passing through the saddle point. 

The contribution to the integral is then composed of three parts 

H = I + I + I. . 1)1 is 2 ui times the residue series 

z p s r> p 

of the poles enclosed between the new and the old paths, 2) 1^ 
is the branch cut contribution due to k^ z and k 2z » and 3) I s 
is the saddle point contribution. Each pole on the 0-plane can 
be interpreted as a normal mode pertaining to the layer medium. 

We therefore discuss our solutions separately for the modal 
part corresponding to the residual series and the non- modal 
part corresponding to the branch cut and the saddle point 
contributions . 

5.2 MODAL SOLUTIONS 

In order to find the normal modes pertaining to the 
layered medium, we determine the poles of the integrand. 

Due to the two branch points at = k^ and = k 2 , there 
are four corresponding Riemann sheets . We are interested 
in the sheet where Re{k lz ) > 0 and Re(k 2z ) > 0. The poles 
of the integrand are in turn determined by the zeros of the 
denominator of 1+R '. Let 

R 12 R 10 ~ e ^ 0 < Re(^) < 2tt (5.4) 

we can write, in view of the denominator of (2) 
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<f> + 2k^ z d = 2 tt i i - 1, 2 , . . . , « 


(5.5) 


For each l, there are two roots corresponding to Re(k z ) > 0 
and Re (kj > 0 . A pole location plot on the complex 6-plane 

t 

is illustrated in Fig. 5.1 for a layer of depth 4A with permi- 
tivity 2 . 16 (1+iO .0022) over a perfect conductor. The 
original path of integration and the new steepest descent path'' 


are also shown on the same figure. 

The poles that lie between the new and the old paths 

of integration contribute to the integral in (5.1). The con- 
tribution is determined by 2iri times the residue series. 


, I£ 

I P T" 


pole: 


r- x 01 (l +Rl2 e i2k l^) 

z 


(5.6) 


d r . , ~ „ i2ki „d. 

hi r tl+E i2 R oi e lz 1 

P 


H { 1) < k p P>e lk * 2 > 


at the pole 


Each term in the summation series corresponds to a normal 
mode. Note that each normal mode is a cylindrical wave. 

Note also that as p <*>, the Hankel function has the asymptotic 

.. .. JT~ ik n p-3iv 

limit ^ P K . 


With regard to the modes that are excited by the antenna, 
we have the following observations: 1) All excited modes 

have Re (k ) < k-. 2) For the class of modes having Re (k ) < k, 

r r 

they all have Re(k ) > 0 and Im(k ) < 0. Physically these 

z z 

modes correspond to waves leaving the surface with increasing 
amplitude. We term these waves leakly modes, they can 
also be called radiation modes of the structure. 3) For the 
class of modes having k < Retk^) < Re(k^) , 


all have 
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Re (k ) < 0 and Im(k ) >0. Physically these nodes correspond 

z z 

to waves coining toward the surface with increasing amplitudes 
away from the surface. We term these waves creeping modes. 

They are surface-like modes and are related to the waveguide 
modes. 4) For the leaky modes, the transverse wave number k p 
has both positive real part arid positive imaginary part. The 
magnitude of the imaginary part is large as compared with 
that of the creeping modes. The leaky modes decay very 
rapidly away from the transmitter and are important only in 
regions near the dipole. 5) The transverse components of 
the creeping modes have a positive real part and a positive 
imaginary part. The creeping waves decay relatively slower 

than the leaky v;aves with trasverse distance. 

The creeping inodes are the significant ones in regions far 
from the transmitting antenna. 6) The number of normal modes 
excited by the antenna depends on the thickness of the slab. 

If the slab is sufficiently thin no mode will be excited 
at all. We discuss the thin layer case in Section 5.4. The calcu- 
lated results are presented in Appendix C for different layer thickness. 

5.3 INTEGRATION AROUND THE SADDLE POINT 


In the absence of any stratified medium, the saddle point 
| contribution from the transmitter corresponds to the direct 

wave from the dipole. When there is stratified medium present, 
the evaluation of the saddle point contribution is complicated 
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by the poles in the neighborhood of the saddle point. Besides 
the pole effects, the saddle point contribution is also affected 
by the rapidly varying exponential exp(i2k^ z d) factor in the 
reflection coefficient. The oscillation of this factor is fast^ 
when d is large. We thus evaluate the saddle point contribution 
for the following three cases: 1) when p >> d, we use the 

ordinary saddle point method (OSP) . 2) When p > d, we use 

the modified saddle point method (MSP) . 3) When pv d, we resort 

to numerical evaluation by using the Gaussian Her mite Guadrature 
(GHQ) formula. 

1) In applying the ordinary saddle point method, we transform 
the integral to the form 

I = (x) exp(“X 2 /2) dx; $(x) = l A ?j x m (5.7) 

m=0 

The solution is then given by 

I s - OSP = /2*{A 0 + A 2 + 3A 4 + . . . ) (5. 8) 

-3 

V?e include terras to p 

The integral (1) can be transformed into the form (5.7) by setting 

. 2 

lkp (1 - sin0) = x /2 (59) 

When p >> d, there are no poles near the saddle point and 
the integrand does not oscillate very fast. Thus the 
ordinary saddle point method is applicable. 

2) Suppose there is a pole near the saddle point at x = Xq. 
Xq = [i2kp(l - sin0p) ] 1//2 


(5.10) 
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where 8^ is the location of the pole. Let 

M = lim $(x) (x - x n ) (5.11) 

x^x 0 o 

VJe can write 

<Mx) » Mx) + (5.12) 

x x 0 


and 


I = W a + W 

S B P 


(5.13) 


where 


W = C — ^ - 

p to ° y.-x. 


‘0 


w o = x q //2 


2 /P 

+iirMe~ x 0' erfcf-iw^) lw 0 >0 


2 n 

-iTTMe"’ x Cr erfc(iWg) lmw^<0 


(5.14) 


(5.15) 


With the effect of the poles being subtracted , W can be 

~ s 


solved asymptotically 


f* 

W = f°° ^( x ) e “ x / 2 dx = OSP - [W J 
S P 


(5.16) 


where (W 3 denotes W in inverse power series of w A . We thus 

p p 0 

have the solution for the modified saddle point method. 

I s « MSP - OSP + po i es (N p - [W p J > ’ (5.17) 

Evidently, as >> 1, the modified saddle point result is 

identical to the ordinary saddle point result. The summation 

in (17) extends over all poles near the saddle point. Both 

-3 

OSP and [V? ] are worked up to p terms. 

3) The Gaussian Hermite quadrature formula is particularly 
useful in evaluating the line integral along the steepest descent 
path. Let 

2 (5.18) 

ikP - ik P = x' 

P 
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we have 


k_ = 


, • 2 

F2 - $ 1 

TTkP 12 IkP J 


( 5 . 19 ) 


And the integral can be cast into the form 


Ig “jr* f(x) exp (-x 2 ) dx 


( 5 . 20 ) 


where 


jri JL * 

4tt k 2 p 


f (X) = Hr- tr— r d + R*") (k p)e iJ, z 


TE X „(1) .. * _ik„z + x 

1 'V 


By the GHQ formula. 


n 


X s - j f 1 w i £<x i> + R „ 


where 


( 5 . 21 ) 


( 5 . 22 ) 


x^ - ith zero of the Hermite polynominl H n (x) 

w. = /tt 2 n ~ 1 nl/n 2 [H , (x . ) ) 2 
l n-1 i 


( 5 . 23 ) 


is the weighting factor, -and 


R - A n! 
n 2 n (2n) 1 


f* 2n) (£) 


-oo < £ < 


( 5 . 24 ) 


is the error term. For the case p <. d, there are poles close 
to the contour of integration , and also the integrand oscillates 
rapidly due to the presence of the exp(i2k^ z d) factor in the 
reflection coefficient. With the knowledge of the locations 
of the poles, we can adjust n such that the error term R is 1 
vanishingly small . 
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Fig. 3 illustrates the case of a 4A dielectric layer 
with permittivity = 2.16(1 + iO .002 ) over a perfect 
reflector. The saddle point contribution is calculated by 
using OSP, MSP and GHQ . It is seen that all three methods 
approach same values as p > 100A. The MSP and GHQ coincides 
for p > 10 A. All three methods depart considerably for 
p < 10 A. Thus for p < 10 A, we recommand GHQ, for 10 A < p < 100 A, 
one can use either MSP or GHQ. For p > 100 A, OSP is sufficient- 
ly accurate. 


5.4 INTEGRATION ALONG THE BRANCHES CUTS 


There are two branch points at k p = k^ and k p = Cor- 

respondingly we choose two branch cuts R s(k lz ) = 0 and 
Re ^ k 2z^ “ The integrand is an even function of k^ z . Thus 

the branch point at k^ does not contribute to the integral. 

The contribution due to the branch point k 2z can be calculated 
by detouring the branch cut to a path parallel to the path of 
the steepest descent. The result is - 

XX k 2 e -2A.2 2 -k 1 2 a B ik 2 P 

p 

T 


(5.25) 


I = H 

■b 2 it 


Ql 10 

i'2k', _d. 


(1+e 1 1 z u Eoi) 


k lz k z 


k = 

P 2 


Physically this corresponds to an inhomogeneous wave that 
decays away from the bottom surface with the exponential factoir 
exp {■ 


(-2v4 2 2 


k l d) • 
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Therefore, except in cases when the slab region is very 
thin, the effect due to the branch point k 2 is not observed 
on the surface of the layered medium. 

In Figure 4, v?e shov? contributions due to the branch point 
for three cases. The permittivity of the first medium is 
3.0(1 + iO.OODeiQ The permittivity of the second medium is 
5.0(i + i0 .002 )£q The depth varies from 0,051, 0.1X, 0.2X, 

0.5 X, 1.0 X. Obviously, as the depth increases, the 1^ effect 
diminishes . 


5.5 RESULTS AND DISCUSSIONS 

In the previous three sections, we have discussed the 
three separate contributions to the integral in (1) . The 
total field at a receiving point due to a transmitting antenna 
is obtained by summing over all the three contributions. With 
applications to the radio frequency interferometry in mind, 
we illustrated the methods by calculating the interference pat- 
terns as a function of distance from the transmitter. 

VThen the slab is sufficiently thin, no mode is excited 
and the integral is due to the saddle point and branch point k 2 « 
Figure 5 shows interference patterns due to layers with 

e. = 3 . 0 (1+iO .00 1) e and depths 0.05X and 0 .IX on top of a 
^ 0 

dielectric half-space with = 5 .0 (1+iC .002} . The saddle point 
contribution is calculated by the ordinary saddle point method. 
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When the layer gets thicker, modes will be excited. We show 
in Figure 6 the interference patterns for the same media with 
depth .2, .5 and 1 wavelength. Only 1 to 3 creeping modes 
are excited. In Figure 7, we illustrated the case of a 4X 
layer on top of a perfect conductor. The pole location plot 
has been shown in Figure 2, Because there are poles that are 
quite close to the saddle point, OSP and MSP are not applicable 
in near and intermediate ranges (Figure 3) . The saddle point 
contribution in Figure 7 has been calculated with GHQ and flSP . 
GHQ result should be more accurate ,e specially in near field 
regions. 

Thus the problem of electromagnetic fields due to a hori- 
zontal electric dipole antenna has been solved with the mode 
approach. We conclude that in calculating the saddle point con- 
tribution in the near field region, one should use the GHQ 
method. In intermediate regions, and in many cases far field 
zones, the MSP method whould be used when there are poles 
near the saddle point. When the layer is thick, all branch point 
contributions can be neglected. But when the layer is very thin, 
one must take into account the branch point contribution. Al- 
though the study is made for the two layer case, the approach 
can be readily extended to more than two layer cases. 



5-12 


Figure 5. 


Figure 5. 


Figure 5. 


Figure 5. 


Figure 5. 


FIGURE CAPTIONS 


1 Pole plot on the complex 0-plane. The depth is 

d — 4Xq , the permittivities are = 2.160 + i0.002. 

The lower half-space is a perfect reflector. 

2 Comparison of results by OSP, MSP and GHQ. The 

depth is d - 4 Aq, the permittivity is = 2 .16 (1+iO . 002) . 
The lower half-space is a perfect reflector. 

3 Branch cut contribution compared for layer depths 0.05 A, 
0.5 A and 1.0A. The permittivities are = 3 .0 (1+iO .001) 
ande^ = 5 .0 (1+iO .002) . The vertical scale is 40db per 
division . 

4 Interference patterns for thin layers . There is no 
pole contributions. The permitivitties are 

e 1 = 3.0 (l+i0.001) and e ^ = 5 .0 (1+iO .002) . The pat- 
terns are plotted for d — 0.05X and d = 0.1 A. Verti- 
cal scale is 20db/division . 

5 Interference patterns for a layer with = 3 .0 (1+iO .001) 
on top of a half-space with - 5 . 0 (1+iO .002) . Calcu- 
lations are done for three layer depths. When d = 0.2A 
and 0.5A, there is only 1 creeping mode being excited. 

Mien d = 1A, 3 creeping modes are excited. Vertical 
scale is 20db per division. 
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Figure 5.6 interference patterns calculated with both MSP and GHQ. 
The depth is d = 4X and the permittivity is 
e «= 2 .6 (1+iO .0022) . The subsurface is a perfect 

conduc tor . 
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Figure 5.3 
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Chapter 6. 
NUMERICAL APPROACH 


6-1 SIMPSON’S RULE 

The integrals in Egs. (1) and (2) can be readily evaluated 
numerically provided they are well defined and properly conver- 
gent. In the actual computation procedures, we change the Hankel 
functions to Bessel functions and integrate from 0 to ® instead 
of from -«> to -Ho. Although the Hankel function has a singularity 
at k p = 0, the integrand as a whole is regular at that point. 

For the integrands to be well defined, we must avoid the branch 
point k - k. When Im(k) ^ 0, the branch point does not lie 

r 

on the path of integration. If Im(k) = 0, we choose the Sommer- 
feld path for the integration. 

Simpson's rule is used in carrying out the numerical inte- 
gration. The amount of computer time needed to obtain a given 
accuracy is reduced in two ways. First, because of the presence 
of the factor exp(ik z z) in the integrand, we choose a non-zero z ; 

the larger the z, the faster is the rate of convergence. Physically, 
this choice corresponds to having the receiving point above the 

surface. In the lunar experiment, the receiving antenna is 
about 2 to 3 meters above ground when mounted on the Lunar 
Roving Vehicle. Secondly, we vary the integration increments 
Ak Q as a function of distance from the transmitter. The magni- 

r ^ 

tude of Akp depends mainly on the rate of oscillation of the 
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Bessel function which for large value of the argument 
proportional to cos k^p. 



is 


For comparison with other techniques, we have evaluated 
numerically the H -component of a two layer model with a 
perfect reflector. The frequency is 8MHz which corresponds to 
a free-space wavelength of Aq= 37.5 meters. The layer has 
dielectric constant 3.3£q, loss tangent 0.01, and a layer thick- 
ness of 4A Q . The results are shown in Fig. 6. The receiving 
antenna has a height of 3 meters. We used Ak p p as small as 
0.5 to determine the increment Ak . With a distance of p<14A_, 

p — 0 

Ak p = 0.5/14Ag ~ 0.004. The computation is stopped when the 
absolute magnitude of the integrand becomes smaller than 0.002 
of the accumulative area. A typical value of the number of 
increments used for Computation is about 5,000. Computation 
time with the IBM 360-65 was about 2 minutes per point as 
shown in Figure 
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6.2 FAST FOURIER TRANSFORM METHOD 

Separate the integral (1) into two parts and I^* 
The first part corresponds to solutions in the absence of 
any stratified medium, and the result is given by the iden- 
tity 


*1 ■ h ■ T- "v - * ikzZ “l 11 < v> ak p 

r —co 2 


( 6 . 1 ) 


The second part includes all effects due to the reflection 
coefficient R^ and after changing to Bessel functions. 


il&k 


- { “ps; 


R TE e lkzZ J, (k p)sin$ 


1' p* 


(6.2) 


which is the integral that we want to solve with the FFT. 

In applying the FFT, we use the formula (Gradshteyn and Ryzhik, 

1965) 


v 


r ®' vkp w ,dk P - ? a - 


/v z + p* 


-) 


(6.3) 


v = v_ + iv_r Re {v + ip} > 0 
R J. 


The integral is written in the following form 


I, - l dk o 9( k o>«"'' Rkp d l< k 0 <» dk 


■2 o P P 


P P 


(6.4) 
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where 


g(k p ) 


i — r TE e ikzZ + V R k £ 

1 4ir k R e K < 

z 


(6.5) 


We can write (Cooley, 1967) for Ak p - 1/2F where F is the 
Nyquist frequency 


I * 1 


9 * k p J NAk 


2 21 


_ _ N 
n 2 


i2xr fk „ 

a(f) e ^for 0 <k < - l)Ak 

P 2 p 


( 6 . 6 ) 


where 


a(f) = / g (k ) e 2lrifk p dk 

—co ^ r 


f = n/NAk 

P 


The factor -2Trf corresponds to in Equation (6) . The right 

hand side of (9) is periodic and does not tend to zero at 

We can multiply (9) by exp(-v R k^) so that the right hand side 

of (9) is sufficiently small for the range of k between 
N 

(^- - 1) Ak and “. Then the final solution becomes, in view 

2 . P 

of the identity (6) 



n=-N/2 


a ( 


n 


NAk 


1 

P 


[1 


(6.7) 

<v r‘sIf 1,/ /V^V' 

P y 
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In choosing the increment Ak p , we recall the two alternatives 

suggested in Section II. If we make the upper half space 

slightly conductive, which corresponds to Im k 0, then 

we choose Ak to be smaller than the distance of the branch 
P 

point from the real axis on the complex k p -plane. If we in- 
sist on a real k, then we choose Ak p such that one of the data 
points coincides with the branch point and such that Ak p is 
smaller than the distance of the pole or branch point closest 
to the real axis on the complex k p -plane. 

To calculate the expansion coefficient a in Equation (9) we 
use the FFT algorithm in a subroutine which, for a given set of 



We first alias the function g(k p ) with period NAk p and denote 
the aliased version by g (k ) . The subroutine requires data 

hr r 

points of g p (k p p) between 0 and (N-l)Ak p . Note that the function 

g (k p) is equal to g(k p ) between 0 and <N/2-l)Ak p but g p (k p ) 

between (|-l)Ak p and (N-l) Ak p is equal to g(k p ) between -(N/2)Ak p 

and -Ak . In the subroutine we calculate 
P 
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N-l 


a ( 


n 


p'NAk J Ak p L t 

p 1=0 




n = 0, 1, . . .N-l 
(6.9) 


and returns the aliased version of a(n/NAk ), a_(n/NAk ). 

P P P 

We must be careful in translating from a^n/UAk^) to 

a(n/NAk ) within the limits ~ 2 X k ~ to 2Ak or< * er to use 

P P 

the result (10) . 

We must choose NAk such that the function g(k ) 

P P 

is sufficiently near zero outside the limits -NAkp/2 to 
NAkp/2. Although mathematically , we could choose z to be 
large, its value is predetermined by the experimental 
arrangement. v R must not be too small. Note that must 
be greater than zero but cannot be larger than z or the ex- 
ponent in (8) will cause g(k p ) to diverge. We choosey = z/2. 

In Figure 1 and Figure 2, we show the interference pat- 
terns for H z calculated with FFT on the IBM 360/65 computer 
for a three layer and a six layer case, respectively. The 

height of the observation point, z = 0.213X which corres- 

0 

ponds to 2 meters for a frequency of 32 MHz (values that were 
used in the lunar experiment) . The computation time was about 
2.5 minutes for both cases. These results are to be compared 
with the results obtained with Simpson's rule and shown also 
in Figures 1 and 2 . 
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6.3 DISCUSSION 

We have shown that with the reflection coefficient 
formulation, the problem of calculating the electromagnetic 
field components due to a dipole transmitting antenna on a 
stratified medium can be solved readily with the FFT algorithm 
and the Simpson's rule and illustrated for for a three 

layer and a six layer case. See Figures 1 and 2. 

In Figure 3, we compare the FFT method with results from 

the geometrical optics approximation and the mode analysis. 

The analytical methods consume very little time, but they are 
applicable only in certain regions according to the nature of 
the approximations involved. The geometrical optics result 
applies when the distance is far from the transmitting antenna 
and when the layer is thick and lossy. The mode method can 
be applied to general cases, and is extremely useful for thin 
layers. The FFT results for some general stratified 
media are also presented in Figures 4-9. 
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FIGURE CAPTIONS 

Figure 1. Comparison of the direct and FFT methods for 

a three layer model in which 

e x = (3.3) (1 + iO. 01) e Q d 1 = IX 

c_ = (5.0) (1 + iO. 02) e d, = 2X 

A 0 <6 

e 3 = (8.0) (1 + i0.04)e 0 

Vertical scale is 10 db per division. 

Figure 2. Comparison of the direct and FFT methods for 
a six layer model in which 


e x = (2.0) (1 + iO . 01) e Q d x = 0.5X 

e 2 = (3.0) (1 + i0.02)e Q d 2 = IX 

e 3 = (4.0) (1 + iO . 03) e Q d 3 = 2X 

= (5.0) (1 + i0.04)e Q d^ — 3X 

e 5 = (6.0) (1 + iO . 05) e Q d 5 = 4X 

e t = (8.0) (1 + iO . 06) e Q 


Vertical scale is 10 db per division. 


Figure 3. Comparison of the FFT , Geometrical Optics (GOA), 
and Mode Theory (MODE) for a two layer model in 
which 

e 1 = 3.3(1 + iO . 02) e Q d 1 = 4X 

e t = 6.0(1 + i0.04)e 0 

Vertical scale is 10 db per division. 
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Chapter 7. 

DISCUSSIONS 

7.1 COMPARISON OF THEORETICAL SOLUTIONS 

In Figure 7.1, we show H as a function of distance cal- 
culated on the basis of the tree different approaches for a sin- 
gle model. The model consists of a single layer, 4A in thickness, 
with e = 33e Q and tanS = 0.01, between free space above and a 
perfect conductor below. The fields are calculated for a receiver 
at a height of 3 meters above surface at 8 MH z . Inspection of 
the figure shows that the results from the mode analysis agree well 
with those obtained from numerical integration. Note in particu- 
lar that the positions of the first peaks in the geometrical optics 
approximation and in the mode analysis occur at different loca- 
tions. In the geometrical optics approach, if this peak is used 
to predict the depth of the subsurface reflector, we see that the 
mode analysis will predict a smaller depth. In the following 
sections, the theoretical calculations are compared with the vari- 
ous experimental results. 

7.2 EXPERIMENTAL CONFIRMATION 

Both the concept of this experiment and the equipment have 
been tested extensively on glaciers and with laboratory-sized 
scale models. Because of the variability of natural materials, 
we used three different glaciers— the Goner Glacier in Switzer- 
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land, the Athabasca Glacier in Alberta, and several glaciers 
that drain the Juneau, Alaska icefields. Each of these glaciers 
had been studied previously. Because the shape, depth, and 
physical properties of each glacier were known already, the data 
obtained on the glaciers can be used to check our theoretical 
expressions. Although we have collected data on several profiles 
on the Gorner, about 50 profiles on the Athabasca and about 120 
profiles on the Juneau glacidrs, with lengths of 1 to 1^ km, 
we shall cite results for only two profiles . These results 
are typical of the other profiles. 

The Gorner glacier, located in southern Switzerland, 
has been studied for more than a century. The rather extensive 
set of data includes the results of gravity, seismic, and (D.C.) 
electrical resistivity surveys. The thickness and shape of the 
glacier were well-determined. Our field gear was quite simple 
and included a General Radio 1330A bridge oscillator for the 
transmitter, a Galaxy R530 communications receiver, and homemade 
antennas. The glacier was quite thick in our test area. In 
Fiq. 7.2 we show the interference pattern of the vertical magnetic 
component for the broadside transmitting antenna at 10 MHz. The 
observed peaks and troughs match very well the corresponding 
features in the theoretical curve calculated for a one layer 
medium with dielectric constant 3.2e Q (l + iO .03) . 

The Athabasca glacier, located about 75 miles south of Jasper, 

Q_1 f) 

Alberta, Canada, has been thoroughly studied also by such other 
methods as seismology, gravity, electrical resistivity measure- 
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ments, and drill holes. The ice thicknesses measured by seismic 
reflections and by direct measurement in boreholes agree quite 
well and we use the profiles reported by Paterson and Savage 
as standards with which to compare the thicknesses determined 
from our electrical sounding technique. In figure 7.3, we show 
typical results, the interference pattern of the vertical magnetic 
component at 2 MHz. Note the excellent match between the experi- 
ment curve and the theortical curves obtained with both the 
geometrical optics approximation and the mode analysis. Our 

interferometry data indicate a depth of 180 meters which is in 
crood agreement with the results of drilling, seismic and gravity 
surveys . 

A scaled model tank, operating at 6 GHz, was used to obtain 

interference patterns for a dielectric layer with dielectric 

constant 2.16e (1 + iO .0022) over an aluminum reflector. The 

o 

general features of the model tank were described by Rossiter 

s 

et al . , and some results from the study of this particular case 
for depths greater than 2A were reported there, with our mode 
analysis, we are now able to match the experimental curves for 
shallower depths. In ficr. 7.4, we show results for layer thick- 
nesses of IX and 2A. The curves are comapred with theoretical 
calculations from both the geometrical optics and the mode ap- 
proaches. Note specifically that for these shallow depths the 
geometrical optics approximation fails to account for even the 
fross features whereas the mode approach fits the experimental 
data excellently. 



7-4 


CONCLUSIONS 

The radiation fields due to a horizontal electric dipole 
laid on the surface of a stratified medium have been calculated 
with three different approaches and compared with the various 
experimental results. The solutions are obtained from the 
ref flection coefficient formulation and written in integral forms. 

In the near field of the transmitting antenna, analytical methods 
involving asymptotic expansion are not applicable. Direct numerical 
integration of the integrals by a computer is the simplest and 
the most useful. The numerical method for near field calculations 
also yields accurate results and uses less computer time than far 
field calculations. When the receiver is far away from the trans- 
mitting antenna, the integrals can be evaluated asymptotically by 
the method of steepest decents. For high lossy media and large 

layer thickness, the geometical optics approach gives rather accurate 
results. The interference patterns calculated from this approach 
can be easily interpreted in terms of ray optics. When losses are 
small and layers are thin, the mode approach is most attractive. 

The results can be interpreted in terms of normal modes of the layered 
medium. Although the calculations and illustrations presented 
in this paper were done for the one layer and the two layer cases, 
the calculations can be readily generalized to handle more layers. 



7-5 


Fig. 7. 


Fig. 7. 


Fig. 7. 


Fig. 7. 


FIGURE CAPTIONS 


1 Comparison of theoretical results. These interference 
patterns are calculated for a dielectric layer with 
dielectric constant 3.3e Q (l + iO .01) overlying a 
perfect conductor. The receiver has a height of 0 .08 
wavelength . 

2 Comparison of Gomer glacier data with theoretical 
results for the vertical magnetic field component H z 
at 10 MHz. The theoretical curve is calculated for 

a one-layer medium with dielectirc constant 3.2e Q {l + i0.03) 

3 A set of Athabasca data taken at 2 MHz, site 3 compared 
with the theoretical results obtained with mode approach 
and geometircal optics approach. The theoretical results 
are calculated for a layer of ice with dielectric constant 
3.3e c ( 1 + iO .15) and depth 1.2X *= 180meters. 

4 Scaled Model Tank experimental data compared with theore- 
tical results obtained with the mode approach and the 
geometrical optics approximation. The model consists of 
a layer of oil with dielectric constant e = 2.16e o 

(1 + i0.0022) . The subsurface reflector is an aluminum 
plate. The upper set of curves is for a layer thickness 
of IX, the lower set for 2X,X being the free space 
wavelength. 
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